如何在 Python 中求解这个微分方程？

Question

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

（+∞ 可以替换为 4）。

---

以上是Falkner-Skan方程。我想得到从 0 到 4 的数值解。

其实我已经看过Numerical Methods in Engineering with Python 3这本书，并且使用了书中的方法，但我无法得到答案。可能书中有一些错误。

这是我的代码：

"""
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

您展示的代码应该是通过将边界值问题简化为初始值问题来实现解决边界值问题的射击方法（https://en.wikipedia.org/wiki/Shooting_method）。 但是，代码中有很多错误。

1. 通过阅读作为解决方案核心的 odeint 函数的文档 (https://docs.scipy.org/doc/scipy-0.18.1/reference/generated/scipy.integrate.odeint.html)，可以发现 odeint 的第一个参数是 func(y, t0, ...)。因此，我们必须将行 def F(x,y): 替换为 def F(y,x):

2. 函数 F(y, x) 应该计算 y 在 x 处的导数。将原始的三阶非线性 ODE 重写为三个一阶 ODE 的系统可以发现： (a) 行 y = np.zeros(3) 没有意义，因为它强制所有导数为零； (b) 它应该是 F1 = y[1] 而不是 F1 = y[2] 并且 (c) 行 F0 = y[0] 不是必需的，可以删除。

3. 我们已经想通了，代码应该实现拍摄方法。这意味着在左边界我们必须设置条件y(0) = 0, y'(0) = 0，这是我们从问题陈述中得到的。但是，要解决初始值问题，我们还需要一个 y''(0) 的条件。我们的想法是，我们将使用 y''(0) = u 形式的不同条件求解我们的 ODE 系统，直到对于 u 的某个值，解决方案在给定容差的右边界处满足边界条件 y'(4) = 1。因此，出于实验目的，让我们将for i in range(1): 行替换为for i in range(5):，并将y'(4) 的值打印为print Y[-1,1]。输出将是：

   -1.26999326625 
   19.263464565
   73.5661968483
   175.047093183
   340.424666137
   

可以看出，如果u = 1 的增量： (a) y'(4) 单调增加和 (b) 在 u = 1 y'(4)=-1.26999326625，在 u = 2 y'(4)=19.263464565。因此，满足边界条件y'(4)=1 的解可以用1<u<2 找到。所以，我们必须减小u的增量，才能找到容差更大的解。

最后，我们可以重写代码如下：

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()

对ustep=0.01 的数值研究表明，有两种不同的解决方案是可能的。一个在0<u<1，另一个在1<u<2。这些解决方案如下所示。

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()