In [23]:
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 notebook  
sns.set()

4th order Runge-Kutta for comparison¶

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

ODE BVP by finite differences¶

In [10]:
import numpy as np
np.set_printoptions(linewidth=300)
In [11]:
def p(x): return 0*x
def q(x): return 100 + 0*x
def phi(x): return 0*x
alpha = 1
beta = 0
x0,x1 = 0,1
In [12]:
def f(x,Y):
    global p,q,phi
    y,yp = Y
    return np.array([ yp, phi(x) -p(x)*yp - q(x)*y ])
In [ ]:
# wild

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)
x0,x1 = 0,1
alpha = 1
beta = 10
sigmaguess = -120
In [ ]:
# sinusoidal y'' + y = 0

def p(x):   return 0*x
def q(x):   return 1 + 0*x
def phi(x): return 0*x 
x0,x1 = 0,3/2*np.pi
alpha = 0
beta = 1
sigmaguess = -1
In [ ]:
# exponential y'' + y' = 0

def p(x):   return 1 + 0*x
def q(x):   return     0*x
def phi(x): return 0*x 
x0,x1 = 0,1
alpha = 2
beta = alpha+1
aa = 1/(1/np.exp(1)-1)
sigmaguess = -aa
In [ ]:
# Golub and Ortega y'' - 100 y = 0   # tough for shooting

def p(x):   return     0*x
def q(x):   return -99 + 0*x
def phi(x): return 1 + 0*x 
x0,x1 = 0,1
alpha = 1
beta  = 0
sigmaguess = -11
In [32]:
# example BVP "wild"

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)
x0,x1 = 0,1   # a,b
alpha = 1
beta = 10

plt.figure(figsize=(8,8))

# compare with RK4
m = 200
h = (x1-x0)/m
for sigmaguess in [-120]:
    ta,Ya = rk4(f,x0,[alpha,sigmaguess],h,m)
    plt.plot(ta,Ya[0,:],label='an accurate RK4-generated solution to compare to')

beta = Ya[0,-1] # change beta to the final condition of the RK4 run so that FD will converge to RK4 trajectory


#N = 5 #99
for N in np.array([5,10,20]):#,40,80,160])-1:
    x = np.linspace(x0,x1,N+2)
    h = x[1]-x[0]
    a = 1/2 - h/4*p(x)
    b = -1 + h**2/2*q(x)
    c = 1/2 + h/4*p(x)
    r = h**2/2*phi(x)
    r[1] -= a[1]*alpha
    r[N] -= c[N]*beta

    # now need to load up the matrix A  - at first as a general matrix 
    A = np.zeros((N,N))
    # let's do it with a loop, though there are slicker ways
    for i in range(1,N+1):
        if i>1: A[i-1,i-2] = a[i]
        A[i-1,i-1] = b[i]
        if i<N: A[i-1,i  ] = c[i]
    #print(A)
    u = np.zeros(N+2)
    u[0] = alpha
    u[-1] = beta
    u[1:-1] = np.linalg.solve(A,r[1:N+1])

    plt.plot(x,u,'.-',label=f'finite differences N={N}')
plt.ylim(-30,30)
plt.legend();
plt.xlabel('x')
plt.ylabel('y')
Out[32]:
Text(0, 0.5, 'y')
In [ ]: