from matplotlib import rcdefaults
rcdefaults() # restore default matplotlib rc parameters
%config InlineBackend.figure_format='retina'
import seaborn as sns # wrapper for matplotlib that provides prettier styles and more
import matplotlib.pyplot as plt # use matplotlib functionality directly
%matplotlib inline
sns.set()
import numpy as np
def rk4(f,t0,y0,h,m):
b21 = 1/2
b31 = 0
b32 = 1/2
b41 = 0
b42 = 0
b43 = 1
b51 = 1/6
b52 = 1/3
b53 = 1/3
b54 = 1/6
t = t0
y = np.array(y0)
ya = np.empty((len(y0),m+1))
ya[:,0] = y
ta = np.linspace(t0,t0+m*h,m+1)
for k in range(m):
t = t0 + k*h
K1 = f(t,y)
y1 = y + h*b21*K1
K2 = f(t+h*b21,y1)
y2 = y + h*(b31*K1 + b32*K2)
K3 = f(t+h*(b31 + b32 ), y2)
y3 = y + h*(b41*K1 + b42*K2 + b43*K3)
K4 = f(t+h*(b41 + b42 + b43 ), y3)
y4 = y + h*(b51*K1 + b52*K2 + b53*K3 + b54*K4)
y = y4
ya[:,k+1] = y
return ta,ya
def f(x,Y):
y,u = Y
return np.array([ u, 3*x*u*y ]) # random nonlinear ODE for illustration
m = 200
y0 = 1 # alpha, the starting value
y1 = 2 # beta, the target ending value
plt.figure(figsize=(12,5))
plt.plot(0,y0,'o')
plt.plot(1,y1,'o')
for u0,person in zip([1],"abcdefg"):
ta,ya = rk4(f,0,[y0,u0],1/m,m)
ya
plt.plot(ta,ya[0,:],alpha = 0.95,label=person+f": y'(0)={u0}")
display( f"{person}: y'(0)={u0}, y(1)={ya[0,-1]}" )
plt.text(1,y1,' target')
plt.xlabel('t',fontsize=15)
plt.ylabel('y(t)\n',fontsize=15,rotation=0); plt.legend(loc='upper center');
"a: y'(0)=1, y(1)=5.133451878500449"