★ Boundary Value Problems ★



In [1]:

    
# Import modules
import numpy as np
import scipy

7.1 Shooting Method

Example

Apply the shooting method to the boundary value problem

$$ \left\{\begin{matrix}\begin{align*} & y'' = 4y \\ & y(0) = 1 \\ & y(1) = 3 \end{align*}\end{matrix}\right. $$



In [383]:

    
def ode_rkf45(f, a, b, y0, h = 1e-3, tol = 1e-6):
    w = y0.astype(np.float64)
    t = a
    while(t < b):
        w_this, t_this = w, t
        s1 = f(t, w)
        hs1 = h * s1
        s2 = f(t + h / 4, w + hs1 / 4)
        hs2 = h * s2
        s3 = f(t + 3 / 8 * h, w + 3 / 32 * hs1 + 9 / 32 * hs2)
        hs3 = h * s3
        s4 = f(t + 12 / 13 * h, w + 1932 / 2197 * hs1 - 7200 / 2197 * hs2 + 7296 / 2197 * hs3)
        hs4 = h * s4
        s5 = f(t + h, w + 439 / 216 * hs1 - 8 * hs2 + 3680 / 513 * hs3 - 845 / 4104 * hs4)
        hs5 = s5 * h
        s6 = f(t + h / 2, w - 8 / 27 * hs1 + 2 * hs2 - 3544 / 2565  * hs3 + 1859 / 4104 * hs4 - 11 / 40 * hs5)
        z = w + h * (16 / 135 * s1 + 6656 / 12825 * s3 + 28561 / 56430 * s4 - 9 / 50 * s5 + 2 / 55 * s6)
        w += h * (25 / 216 * s1 + 1408 / 2565 * s3 + 2197 / 4104 * s4 - s5 / 5)
        t += h
        e = abs(z - w)
        if np.all(e / np.abs(w) < tol):
            w = z
        else:
            h = 0.8 * pow(tol * abs(w) / e, 1 / 5) * h
        
    return w



In [384]:

    
from scipy.optimize import fsolve

f = lambda t, w : np.array([w[1], 4 * w[0]])

def F(s):
    return (ode_rkf45(f, 0, 1, np.array([1, s])) - 3)[0]

x = fsolve(F, 0)
print(x)









    



[-0.42030605]



In [385]:

    
from scipy.integrate import ode
from scipy.optimize import fsolve

def F(s):
    a = 0
    ya = 1.0
    b = 1
    step = 100
    dt = (b - a) / step
    f = lambda t, w : np.array([w[1], 4 * w[0]])
    r = ode(f).set_integrator('dopri5').set_initial_value(np.array([ya, s]), 0)
    while r.successful() and step > 0:
        x = r.integrate(r.t + dt)[0]
        step -= 1
        
    return x - 3

x = fsolve(F, 0)
print(x)









    



[-0.42030605]

7.2 Finite Difference Methods

$$ \begin{align*} y'(t) &= \frac{y(t+h) - y(t-h)}{2h} - \frac{h^2}{6}y'''(c) \\ y''(t) &= \frac{y(t+h) - 2y(t) + y(t-h)}{h^2} + \frac{h^2}{12}f''''(c) \end{align*} $$

Example

Solve the BVP

$$ \left\{\begin{matrix}\begin{align*} y'' &= 4y \\ y(0) &= 1 \\ y(1) &= 3 \end{align*}\end{matrix}\right. $$

using finite differences

$$ \begin{align*} & \frac{w_{i+1} - 2w_{i} + w_{i-1}}{h^2} - 4w_i = 0 \\ & \rightarrow w_{i-1} + (-4h^2 - 2)w_i + w_{i_+1} = 0 \\ & \Rightarrow \\ & 1 + (-4h^2 - 2)w_1 + w_2 = 0 \\ & w_1 + (-4h^2 - 2)w_2 + w_3 = 0 \\ & w_2 + (-4h^2 - 2)w_3 + 3 = 0 \\ \end{align*} $$$$ \begin{bmatrix} -\frac{9}{4} & 1 & 0 \\ 1 & -\frac{9}{4} & 1 \\ 0 & 1 & -\frac{9}{4} \end{bmatrix} \begin{bmatrix} w_1 \\ w_2 \\ w_3 \end{bmatrix} = \begin{bmatrix} -1 \\ 0 \\ -3 \end{bmatrix} $$



In [5]:

    
A = np.array([-9/4, 1, 0, 1, -9/4, 1, 0, 1, -9/4]).reshape(3, 3)
b = np.array([-1, 0, -3])
x = np.linalg.solve(A, b)
print(x)









    



[ 1.02494331  1.30612245  1.9138322 ]

7.3 Collocation and the Finite Element Method

Example

Apply the Finite Element Method to the BVP

$$ \left\{\begin{matrix}\begin{align*} y'' &= 4y \\ y(0) &= 1 \\ y(1) &= 3 \end{align*}\end{matrix}\right. $$



In [50]:

    
def finite_element_method_sp( a, b, ya, yb, n):
    h = (b - a) / (n + 1)
    alpha = (8 / 3) * h + 2 / h
    beta = (2 / 3) * h - 1 / h
    A = np.zeros(n * n).reshape(n, n)
    np.fill_diagonal(A, alpha)
    dia_range = np.arange(n - 1)
    A[dia_range, dia_range + 1] = beta
    A[dia_range + 1, dia_range] = beta
    b = np.zeros(n)
    b[0] = -ya * beta
    b[-1] = -yb * beta
    x = np.linalg.solve(A, b)
    return x
    
finite_element_method_sp(0, 1, 1, 3, 3)









    Out[50]:





array([ 1.01091223,  1.2855407 ,  1.89552762])