How to solve this differential equation in Python?

Question

y''' + yy" + (1 - y'^2)= 0 y(0)=0, y'(0)=0, y'(+∞)=1

(The +∞ can be replaced with 4).

---

The above is a Falkner-Skan equation. I want to get the numerical solution from 0 to 4.

Actually, I have read the book Numerical Methods in Engineering with Python 3, and used the methods in the book, but I cannot get the answer. Perhaps there are some errors in the book.

Here is my code:

"""
solve the differential equation
y''' + yy" + (1 - y'^2)= 0  y(0) = 0 y'(0) = 0,y'(+inf) = 0
"""
import numpy as np
from scipy.integrate import odeint
from printSoln import *
start = time.clock()
def F(x,y):   # First-order differential equations
    y = np.zeros(3)
    F0 = y[0]
    F1 = y[2]
    F2 = y[2]
    F3 = -y[0]*y[2] - (1 - y[1]*y[1])
    return [F1,F2,F3]

def init(u):
    return np.array([0,0,u])
X = np.linspace(0, 5, 50)
u = 0
for i in range(1):
    u += 1
    Y = odeint(F, init(u), X)
    print(Y)
    if abs(Y[len(Y)-1,1] - 1) < 0.001:
        print(i,'times iterations')
        printSoln(X,Y,freq=2)
        print("y'(inf)) = ", Y[len(Y)-1,1])
        print('y"(0) =',u)
        break
end = time.clock()

print('duration =',end - start)

Answer 1

The code you show is supposed to realize the shooting method to solve boundary value problem by reducing it to the initial value problem (https://en.wikipedia.org/wiki/Shooting_method). However, there are many errors in the code.

1. By reading documentation (https://docs.scipy.org/doc/scipy-0.18.1/reference/generated/scipy.integrate.odeint.html) for odeint function which is a core of the solution, one can find that first argument of the odeint is a func(y, t0, ...). Therefore, we have to replace the line def F(x,y): by def F(y,x):

2. Function F(y, x) should computes the derivative of y at x. Rewriting the original 3rd order non-linear ODE as a system of three 1st order ODEs one can find that: (a) line y = np.zeros(3) does not make sense, since it forces all derivatives to be zero; (b) it should be F1 = y[1] instead of F1 = y[2] and (c) line F0 = y[0] is not necessary and can be deleted.

3. We have already figured out that the code should realize shooting method. It means that at the left boundary we have to set conditions y(0) = 0, y'(0) = 0 which we have from the problem statement. However, to solve initial value problem we need one more condition for y''(0). The idea is that we will solve our system of ODEs with different conditions of the form y''(0) = u until for some value of u the solution satisfies boundary condition y'(4) = 1 at the right boundary with given tolerance. Hence, for experiment purposes, let's replace line for i in range(1): by for i in range(5): and print the value of y'(4) by print Y[-1,1]. The output will be:

   -1.26999326625 
   19.263464565
   73.5661968483
   175.047093183
   340.424666137
   

It can be seen that if increment of u = 1: (a) y'(4) monotonically increase and (b) at u = 1 y'(4)=-1.26999326625, at u = 2 y'(4)=19.263464565. Therefore, the solution which satisfies boundary condition y'(4)=1 can be found with 1<u<2. So, we have to decrease the increment of u in order to find the solution with higher tolerance.

Finaly, we can rewrite the code as follows:

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint
import time

start = time.clock()
def F(y, x):   # First-order differential equations
    F1 = y[1]
    F2 = y[2]
    F3 = -y[0]*y[2] - (1 - y[1]*y[1])
    return [F1,F2,F3]

def init(u):
    return np.array([0,0,u])

X = np.linspace(0, 4, 50)
niter = 100
u = 0
tol = 0.1
ustep = 0.01

for i in range(niter):
    u += ustep
    Y = odeint(F, init(u), X)
    print Y[-1,1]
    if abs(Y[len(Y)-1,1] - 1) < tol:
        print(i,'times iterations')
        print("y'(inf)) = ", Y[len(Y)-1,1])
        print('y"(0) =',u)
        break

end = time.clock()
print('duration =',end - start)

plt.plot(X, Y[:,0], label='y')
plt.plot(X, Y[:,1], label='y\'')
plt.legend()
plt.show()

Numerical investigations with ustep=0.01 shows that two diferent solutions are possible. One at 0<u<1 and another at 1<u<2. These solutions are shown below.

Answer 2

# biscetion
""" root = bisection(f,x1,x2,switch=0,tol=1.0e-9).
Finds a root of f(x) = 0 by bisection.
The root must be bracketed in (x1,x2).
Setting switch = 1 returns root = None if
f(x) increases upon bisection.
"""
import math
from numpy import sign
def bisection(f,x1,x2,switch=1,tol=1.0e-9):
    f1 = f(x1)
    if f1 == 0.0: return x1
    f2 = f(x2)
    if f2 == 0.0: return x2
    if sign(f1) == sign(f2):
        print('Root is not bracketed')
    n = int(math.ceil(math.log(abs(x2 - x1)/tol)/math.log(2.0)))
    for i in range(n):
        x3 = 0.5*(x1 + x2); f3 = f(x3)
        if (switch == 1) and (abs(f3) > abs(f1)) \
            and (abs(f3) > abs(f2)):
            return None

        if f3 == 0.0: return x3
        if sign(f2)!= sign(f3): x1 = x3; f1 = f3
        else: x2 = x3; f2 = f3
    return (x1 + x2)/2.0
""""//////////////////////////////////////////////////////////////////"""
# run_kut4
""" X,Y = integrate(F,x,y,xStop,h).
    4th-order Runge-Kutta method for solving the
    initail value problem {y}' = {F(x,{y})},where
    {y} = {y[0],y[1],...,y[n-1]}.
    x,y = initial conditions
    xStop = terminal value of x
    h     =increment of x used in integration
    F     =user - supplied function that returns the 
           array F(x,y) = {y'[0],y'[1],...,y'[n-1]}.
"""
import numpy as np
def integrate(F,x,y,xStop,h):

    def run_kut4(F,x,y,h):
        K0 = h*F(x,y)
        K1 = h*F(x + h/2.0, y + K0/2.0)
        K2 = h*F(x + h/2.0, y + K1/2.0)
        K3 = h*F(x + h, y + K2)
        return (K0 + 2.0*K1 + 2.0*K2 + K3)/6.0
    X = []
    Y = []
    X.append(x)
    Y.append(y)
    while x < xStop:
        h = min(h,xStop - x)
        y = y + run_kut4(F,x,y,h)
        x = x + h
        X.append(x)
        Y.append(y)
    return np.array(X), np.array(Y)
""""//////////////////////////////////////////////////////////////////"""    
## printsSoln
""" printSoln(X,Y,freq).
    Prints X and Y returner from the differential
    equation solves using printout frequency 'freq'.
        freq = n prints every nth step.
        freq = 0 prints inital and final values only.
"""
def printSoln(X,Y,freq):

    def printHead(n):
        print('\n{:>13}'.format('x'),end=" ")
        for i in range(n):
            print('{}{}{}{}'.format(' '*(9),'y[',i,']'),end=" ")
        print()
    def printLine(x,y,n):
        print("{:13.4}".format(x),end=" ")
        for i in range(n):
            print("{:13.4f}".format(y[i]),end=" ")
        print()
    m = len(Y)
    try: n  = len(Y[0])
    except TypeError: n = 1
    if freq == 0: freq = m
    printHead(n)
    for i in range(0,m,freq):
        printLine(X[i],Y[i],n)
    if i != m - 1: printLine(X[m - 1],Y[m - 1],n)
""""//////////////////////////////////////////////////////////////////""" 

"""
solve the differential equation
y''' + yy" + B(1 - y'^2)= 0  y(0) = 0 y'(0) = 0,y'(+inf) = 0
B = [0,0.2,0.4,0.6,0.8,1.0,-0.1,-0.15]
"""
import matplotlib.pyplot as plt
import time
start = time.clock()
def initCond(u): #Init. values of [y.y']; use 'u' if unkown 
    return np.array([0.0,0.0,u])

def r(u):   # Boundart condition resdiual
    X,Y = integrate(F,xStart,initCond(u),xStop,h)
    y = Y[len(Y) - 1]
    r = y[1] - 1.0
    return r

def F(x,y):   # Third-order differential equations
    F = np.zeros(3)
    F[0] = y[1]
    F[1] = y[2]
    F[2] = -y[0]*y[2] - (1 - y[1]*y[1])
    return F

xStart = 0.0   # Start of integration
xStop = 5    # End of integraiton
u1 = 1.0       # 1st trial value of unknown init. cond.
u2 = 2.0       # 2nd trial value of unknown init. cond.
h = 0.1        # Step size
freq = 2       # Printout frequency


u = bisection(r,u1,u2,tol = 1.0e-9) #Comput the correct initial condition
X,Y = integrate(F,xStart,initCond(u),xStop,h)
print(initCond(u))
printSoln(X,Y,freq)
end = time.clock()
print('duration =',end - start)
plt.plot(X, Y[:,0], label='y')
plt.plot(X, Y[:,1], label='y\'')
plt.legend()
plt.show()