Shooting method for BVP - Numerical Analysis

def f(u):
    rhs = zeros(4)
    rhs[0] = u[1]
    rhs[1] = -u[0] + 2.0*u[1]**2/u[0]
    rhs[2] = u[3]
    rhs[3] = (-1.0-2.0*(u[1]/u[0])**2)*u[2] + 4.0*u[1]*u[3]/u[0]
    return rhs

def initialcondition(s):
    u = zeros(4)
    u[0] = 1.0/(exp(1)+exp(-1))
    u[1] = s
    u[2] = 0.0
    u[3] = 1.0
    return u

def yexact(x):
    return 1.0/(exp(x) + exp(-x))

The function below implements the two step RK scheme.

def solvephi(n,s):
    h = 2.0/n
    u = zeros((n+1,4))
    u[0] = initialcondition(s)
    for i in range(1,n+1):
        u1 = u[i-1] + 0.5*h*f(u[i-1])
        u[i] = u[i-1] + h*f(u1)
    phi = u[n][0] - 1.0/(exp(1)+exp(-1))
    dphi= u[n][2]
    return phi,dphi,u

Let us see the solution corresponding to the initial guess.

n = 100
s = 0.2
p, dp, u = solvephi(n,s)
x = linspace(-1.0,1.0,n+1)
xe = linspace(-1.0,1.0,1000)
ye = yexact(xe)
plot(x,u[:,0],xe,ye)
grid(True), title("s = " + str(s))
legend(("Numerical","Exact"));

The initial guess is far from the solution. We will not start the shooting method.

def shoot(s,n,plotsol=False):
    maxiter, TOL = 100, 1.0e-8
    x  = linspace(-1.0,1.0,n+1)
    if plotsol:
        figure(figsize=(8,5))
        xe = linspace(-1.0,1.0,1000)
        ye = yexact(xe)
        plot(xe,ye,label="Exact")
    it = 0
    while it < maxiter:
        p, dp, u = solvephi(n,s)
        if plotsol:
            plot(x,u[:,0],label=str(it))
        if abs(p) < TOL:
            break
        s = s - p/dp
        it += 1
        print('%5d %16.6e %16.6e'%(it,s,abs(p)))
    if plotsol:
        legend(); grid(True)
    return s,p,x,u

Solve with initial guess of $s=0.2$ and step size in RK method of $h = 2/n$ where $n=64$ .

s, n = 0.2, 64
s,p,x,u = shoot(s,n,True)
print("Root s        = %e" % s)
print("phi(s)        = %e" % p)
sexact = (exp(1) - exp(-1))/(exp(1) + exp(-1))**2
print("Error in root = ",  sexact-s)
figure(figsize=(8,5))
plot(x,u[:,0],xe,ye)
grid(True)
title("Final solution, s = " + str(s))
legend(("Numerical","Exact"));

    1     2.712231e-01     1.113574e-01
    2     2.534535e-01     1.224071e-01
    3     2.472480e-01     2.633274e-02
    4     2.467467e-01     1.840435e-03
    5     2.467438e-01     1.033433e-05
Root s        = 2.467438e-01
phi(s)        = 3.295230e-10
Error in root =  3.335992702416246e-05

Numerical Analysis

ODE with periodic solution

Numerical Analysis

Nonlinear BVP using finite difference method