布拉休斯方程數值求解

布拉休斯方程如下：

\[\begin{equation} f f^{''}+2f^{'''}=0 \\ f(0)=f^{'}(0)=0;f^{''}(0)=1 \end{equation} \]

這是一個非線性常微分方程，下面我們利用四階龍格庫塔方法來求解該方程。

我們引入新的變量：

\[\begin{equation} y_1=f \\ y_2=f^{'} \\ y_3=f^{''} \end{equation} \]

則布拉休斯方程可以等價的表示乘如下：

\[\begin{equation} \left\{ \begin{aligned} &y_1^{'} = y_2 \\ &y_2^{'} = y_3 \\ &y_3^{'} = -\frac{1}{2}y_1 y_3 \end{aligned} \right . \end{equation} \]

利用四階龍哥庫塔有：

\[\begin{equation} \left\{ \begin{aligned} &y_{n+1}=y_n+\frac{h}{6}(K_1+2K_2+2K_3+K_4) \\ &K_1 = f(x_n, y_n) \\ &K_2 = f(x_n+\frac{h}{2}, y_n+\frac{h}{2}K_1) \\ &K_3 = f(x_n+\frac{h}{2}, y_n+\frac{h}{2}K_2) \\ &K_4 = f(x_n+h, y_n+h K_3) \\ \end{aligned} \right . \end{equation} \]

對於（3）式中的第一式利用龍哥庫塔方法，有：

\[\begin{equation} \left\{ \begin{aligned} &y_1^{n+1}=y_1^{n}+\frac{h}{6}(K_{11}+2K_{12}+2K_{13}+K_{14}) \\ &K_{11} = y_2^{n} \\ &K_{12} = y_2^{n} + \frac{h}{2}K_{11} \\ &K_{13} = y_2^{n} + \frac{h}{2}K_{12} \\ &K_{14} = y_2^{n} + h K_{13} \\ \end{aligned} \right . \end{equation} \]

同理，對於(3)中的第二式：

\[\begin{equation} \left\{ \begin{aligned} &y_2^{n+1}=y_2^{n}+\frac{h}{6}(K_{21}+2K_{22}+2K_{23}+K_{24}) \\ &K_{21} = y_3^{n} \\ &K_{22} = y_3^{n} + \frac{h}{2}K_{21} \\ &K_{23} = y_3^{n} + \frac{h}{2}K_{22} \\ &K_{24} = y_3^{n} + h K_{23} \\ \end{aligned} \right . \end{equation} \]

對於(3)中的第三式：

\[\begin{equation} \left\{ \begin{aligned} &y_3^{n+1}=y_3^{n}+\frac{h}{6}(K_{31}+2K_{32}+2K_{33}+K_{34}) \\ &K_{31} = -\frac{1}{2}y_1^{n}y_3^{n} \\ &K_{32} = -\frac{1}{2}(y_1^{n} + \frac{h}{2}K_{31})(y_3^{n} + \frac{h}{2}K_{31}) \\ &K_{33} = -\frac{1}{2}(y_1^{n} + \frac{h}{2}K_{32})(y_3^{n} + \frac{h}{2}K_{32}) \\ &K_{34} = -\frac{1}{2}(y_1^{n} + hK_{33})(y_3^{n} + hK_{33}) \\ \end{aligned} \right . \end{equation} \]

編寫代碼：

import numpy as np

# 四階龍格庫塔（RK4）
L = 10     # length of boundary
h = 0.1    # step length
N = int(L / h)
y1 = np.zeros(N)
y2 = np.zeros(N)
y3 = np.zeros(N)
y3[0] = 0.33206   # test value

for ii in range(N-1):
    k11 = y2[ii]
    k21 = y3[ii]
    k31 = -1 / 2 * y1[ii] * y3[ii]
    k12 = y2[ii] + h / 2 * k21
    k22 = y3[ii] + h / 2 * k31
    k32 = -1 / 2 * (y1[ii] + h / 2 * k11) * (y3[ii] + h / 2 * k31)
    k13 = y2[ii] + h / 2 * k22
    k23 = y3[ii] + h / 2 * k32
    k33 = -1 / 2 * (y1[ii] + h / 2 * k12) * (y3[ii] + h / 2 * k32)
    k14 = y2[ii] + h * k22
    k24 = y3[ii] + h * k32
    k34 = -1 / 2 * (y1[ii] + h * k13) * (y3[ii] + h * k33)
    y1[ii+1] = y1[ii] + h / 6 * (k11 + 2 * k12 + 2 * k13 + k14)
    y2[ii+1] = y2[ii] + h / 6 * (k21 + 2 * k22 + 2 * k23 + k24)
    y3[ii+1] = y3[ii] + h / 6 * (k31 + 2 * k32 + 2 * k33 + k34)

利用上面計算的結果，繪制圖像如下：

import matplotlib.pyplot as plt
x = np.arange(N) * h
fig = plt.figure(figsize=(16, 9))
plt.plot(x, y2, '--', color='red')
# plt.xlabel('$\eta=(u_e/\nu x)^0.5$')
# plt.ylabel("$f^'=u/u_e$")
plt.xlabel('x')
plt.ylabel('f')
plt.grid()
plt.show()

圖像如下：

在上面計算中，我們給出一個嘗試值y3[0]=0.33206。但是，其實我們並不知道\(f^{''}(0)\)的值。因此采用打靶的思想進行確定。具體思路是給出一個\(f^{''}(0)\)的值，計算\(f\)的值，根據插值來修正。（我們知道精確值為\(f(\infty)=1\)）。具體代碼如下：

# 采用修正算法來自動尋找初始值。

def RK4(ic):
    # 四階龍格庫塔（RK4）
    L = 40      # length of boundary
    h = 0.1     # step length
    N = int(L / h)
    y1 = np.zeros(N)
    y2 = np.zeros(N)
    y3 = np.zeros(N)
    y3[0] = ic

    for ii in range(N-1):
        k11 = y2[ii]
        k21 = y3[ii]
        k31 = -1 / 2 * y1[ii] * y3[ii]
        k12 = y2[ii] + h / 2 * k21
        k22 = y3[ii] + h / 2 * k31
        k32 = -1 / 2 * (y1[ii] + h / 2 * k11) * (y3[ii] + h / 2 * k31)
        k13 = y2[ii] + h / 2 * k22
        k23 = y3[ii] + h / 2 * k32
        k33 = -1 / 2 * (y1[ii] + h / 2 * k12) * (y3[ii] + h / 2 * k32)
        k14 = y2[ii] + h * k22
        k24 = y3[ii] + h * k32
        k34 = -1 / 2 * (y1[ii] + h * k13) * (y3[ii] + h * k33)
        y1[ii+1] = y1[ii] + h / 6 * (k11 + 2 * k12 + 2 * k13 + k14)
        y2[ii+1] = y2[ii] + h / 6 * (k21 + 2 * k22 + 2 * k23 + k24)
        y3[ii+1] = y3[ii] + h / 6 * (k31 + 2 * k32 + 2 * k33 + k34)
    return y2[-1]

epsilon = 1e-9
# 開始迭代測試修正
ic = 1
alpha = 0.3
for ii in range(20):
    print(ic)
    y3 = RK4(ic)
    error = y3 - 1.0
    if np.abs(error) < epsilon:
        break
    else:
        ic -= alpha * error

運行上面的程序，我們很容易得到\(f^{''}(0)\)從初始值（嘗試值1）降低到正確值0.33205。

1
0.6743617428806606
0.4932473770697171
0.4026815319928927
0.36152169076308643
0.34402498259791453
0.3368568758014264
0.3339705495653165
0.33281689350656146
0.33235717556926003
0.3321742068752738
0.3321014204467237
0.33207247103980964
0.3320609578581748
0.3320563792056564
0.33205455835342595
0.33205383423498663
0.33205354626735756
0.33205343174839597
0.3320533862065114