BVP by shooting¶

imports¶

In [3]:
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 [4]:
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 [5]:
def f(x,Y):
    global p,q,phi
    y,u = Y
    return np.array([ u, phi(x) -p(x)*u - q(x)*y ])


# coefficiencts of some random example 2nd order linear ODE
def p(x):   return 5*x**2
def q(x):   return 30+50*np.sin(30*x)
def phi(x): return 50*np.cosh(x)

Shooting from t=0 to t=1¶

In [12]:
m = 200
y0 = 18  # alpha, the starting value
y1 = -20  # 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, -100 ],["Leeyah's choice","Yifan's choice"]):
    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)',fontsize=15,rotation=0); plt.legend(loc='upper center');

"Leeyah's choice: y'(0)=1, y(1)=8.384923184049143"
"Yifan's choice: y'(0)=-100, y(1)=15.27443731049142"

Can we now come up with the right value of y'(0)?¶

In [7]:
rhostar = (y1-8.384923184049143)/(15.27443731049142-(8.384923184049143))
rhostar
u0 = (1-rhostar)*1 + rhostar*(-100)
u0
Out[7]:
417.1218322473212
In [10]:
m = 200
y0 = 18  # alpha, the starting value
y1 = -20  # 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, -100, 417.1218322473212 ],["Leeyah's choice","Yifan's choice","Calculated 3rd choice"]):
    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)',fontsize=15,rotation=0); plt.legend();

"Leeyah's choice: y'(0)=1, y(1)=8.384923184049143"
"Yifan's choice: y'(0)=-100, y(1)=15.27443731049142"
"Calculated 3rd choice: y'(0)=417.1218322473212, y(1)=-19.999999999999964"
In [14]:
from mpl_toolkits.mplot3d import Axes3D
%matplotlib notebook
fig = plt.figure(figsize=(10,7))
ax = fig.gca(projection='3d')
n=20
ya3 = np.empty((2,m+1,n))
for k,yp0 in enumerate( np.linspace(-30,30,n) ):
    ta,ya = rk4(f,0,[y0,yp0],1/m,m)
    ax.plot(ta,ya[0],ya[1],alpha = 0.95)
    ya3[:,:,k] = ya
ax.plot([0]*n,ya3[0, 0,:],ya3[1, 0,:],color='g')
ax.plot([1]*n,ya3[0,-1,:],ya3[1,-1,:],color='r')
ax.set_xlabel('t')
ax.set_ylabel('y')
ax.set_zlabel('y \'')

/tmp/ipykernel_51163/543498099.py:4: MatplotlibDeprecationWarning: Calling gca() with keyword arguments was deprecated in Matplotlib 3.4. Starting two minor releases later, gca() will take no keyword arguments. The gca() function should only be used to get the current axes, or if no axes exist, create new axes with default keyword arguments. To create a new axes with non-default arguments, use plt.axes() or plt.subplot().
  ax = fig.gca(projection='3d')
Out[14]:
Text(0.5, 0, "y '")
In [ ]:
In [ ]: