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()


    

    















    












    




/home/adam/.local/lib/python3.6/site-packages/ipykernel_launcher.py:9: RuntimeWarning: invalid value encountered in sqrt
  if __name__ == '__main__':
/home/adam/.local/lib/python3.6/site-packages/ipykernel_launcher.py:10: RuntimeWarning: invalid value encountered in sqrt
  # Remove the CWD from sys.path while we load stuff.






    















    

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 [ ]: