布拉休斯方程数值求解

布拉休斯方程如下:

\[\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
posted @ 2020-04-03 09:04  ixum  阅读(578)  评论(0编辑  收藏  举报