In [1]:
%run Steppers.ipynb #Import stepper functions (Euler, RK4, etc)
%run ODEs.ipynb #Import ODEs (Pendulum etc)
%run Tools.ipynb
import matplotlib
matplotlib.use('nbagg')
import matplotlib.pyplot as plt
In [2]:
def SolarDynamo(f,t,a,b):
x,y=f
dxdt=-a*x+2*x*y
dydt=b-x**2-y**2
return(np.array([dxdt,dydt]))
In [3]:
def XNullclines(x,a,b):
y=x*0+a/2
return y
def XNullclineZero(y,a,b):
return(y*0)
def YNullclines(x,a,b):
y1=np.sqrt(b-x**2)
y2=-np.sqrt(b-x**2)
return y1,y2
In [4]:
t0 = 0
tf = 50
h = 0.005
t = np.arange(t0,tf,h) #Generate timesteps
a=0
b=1
xLims,yLims=(-1.2*np.sqrt(b),1.2*np.sqrt(b)),(-1.2*np.sqrt(b),1.2*np.sqrt(b))
x0=0.5
y0=0.5
f0=(x0,y0)
sol=RK4(SolarDynamo, f0, t,args=(a,b))
x,y=np.hsplit(sol,len(sol[0]))
plt.title(r"Solar Dynamo system $a={}, b={}$".format(a,b))
plt.plot(t,x,label="x(t)")
plt.plot(t,y,label="y(t)")
plt.xlabel("Time (s)")
plt.legend()
plt.show()
for f0 in [(0.5,0.5),(0.5,-0.5),(-0.5,0.5),(-0.5,-0.5)]:
sol=RK4(SolarDynamo, f0, t,args=(a,b))
x,y=np.hsplit(sol,len(sol[0]))
plt.plot(x,y, color="C0")
xPlot=np.linspace(-1.05*np.sqrt(b),1.05*np.sqrt(b),5000)
yPlot=np.linspace(-1.05*np.sqrt(b),1.05*np.sqrt(b),5000)
y1=XNullclines(xPlot,a,b)
y2,y3=YNullclines(xPlot,a,b)
x1=XNullclineZero(yPlot,a,b)
plt.plot(xPlot,y1,label="X Nullclines",color="C1")
plt.plot(xPlot,y2,label="Y Nullclines",color="C2")
plt.plot(xPlot,y3,color="C2")
plt.plot(x1,yPlot,color="C1")
plt.xlim(xLims)
plt.ylim(yLims)
PlotQuiver(SolarDynamo,xLims,yLims, args=(a,b))
plt.xlabel("x")
plt.ylabel("y")
plt.legend()
plt.title("Phase Space Plot")
plt.show()
In [5]:
b=np.linspace(0,5,1000)
a=4*np.sqrt(b)
plt.plot(a,b, label=r"$\tau^2-4\Delta=0$")
a=2*np.sqrt(b)
plt.plot(a,b,label=r"$\Delta=0$")
plt.xlabel("a")
plt.ylabel("b")
plt.show()
b=np.linspace(0,5,1000)
a=4*np.sqrt(b)
plt.plot(a,b, label=r"$\tau^2-4\Delta=0$")
a=-2*np.sqrt(b)
plt.plot(a,b,label=r"$\Delta=0$")
plt.xlabel("a")
plt.ylabel("b")
plt.show()
In [ ]: