A nonlinear BVP by shooting¶

imports¶

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

4-stage 4th order Runge-Kutta¶

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

A 2nd order scalar DE expressed as 2D first order system¶

In [19]:
def f(x,Y):
    y,u = Y
    return np.array([ u, 3*x*u*y ])  # random nonlinear ODE for illustration

Shooting from t=0 to t=1¶

In [35]:
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"

Now let's solve using Newton's method, solving the variational equation in parallel to get "$g'$"¶

In [ ]: