Consider the BVP
y′′=−y+y2(y′)2,−1<x<+1 with boundary conditions
y(−1)=y(+1)=(e+e−1)−1 The exact solution is
y(x)=(ex+e−x)−1 Formulate as an initial value problem
y′′=−y+y2(y′)2,−1<x<+1 with initial conditions
y(−1)=(e+e−1)−1,y′(−1)=s Let the solution of this IVP be denoted as y(x;s). We have to determine s as the root of the function
ϕ(s)=y(1;s)−(e+e−1)−1=0 The exact root is
s=y′(−1)=(e+e−1)2e−e−1≈0.24677717 Newton method
To find the root of ϕ(s), we use Newton method with an initial guess s0. Then the Newton method updates the guess by
sm+1=sm−dsdϕ(sm)ϕ(sm),m=0,1,2,… Define
zs(x)=∂s∂y(x;s) Then
zs′(x)=∂s∂y′(x;s) The derivative of ϕ(s) is given by
dsdϕ(s)=∂s∂y(1;s)=zs(1) We can write an equation for zs(x)
zs′′=[−1−2(yy′)2]zs+4yy′zs′ with initial conditions
zs(−1)=0,zs′(−1)=1 First order ODE system
Define the vector
u=⎣⎡u1u2u3u4⎦⎤=⎣⎡yy′zszs′⎦⎤ Then
u1′u2′u3′u4′=u2=−u1+u12u22=u4=[−1−2(u1u2)2]u3+4u1u2u4 Hence we get the first order ODE system
u′=f(u)=⎣⎡u2−u1+u12u22u4[−1−2(u1u2)2]u3+4u1u2u4⎦⎤ with initial condition
u(−1)=⎣⎡(e+e−1)−1s01⎦⎤ Once we solve this IVP, we get
ϕ(s)=u1(1)−(e+e−1)−1,dsdϕ(s)=zs(1)=u3(1) We now code the problem specific data.
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
