布拉休斯方程如下：
[egin{equation}
f f^{''}+2f^{'''}=0 \
f(0)=f^{'}(0)=0;f^{''}(0)=1
end{equation}
]
这是一个非线性常微分方程，下面我们利用四阶龙格库塔方法来求解该方程。
我们引入新的变量：
[egin{equation}
y_1=f \
y_2=f^{'} \
y_3=f^{''}
end{equation}
]
则布拉休斯方程可以等价的表示乘如下：
[egin{equation}
left{
egin{aligned}
&y_1^{'} = y_2 \
&y_2^{'} = y_3 \
&y_3^{'} = -frac{1}{2}y_1 y_3
end{aligned}
ight .
end{equation}
]
利用四阶龙哥库塔有：
[egin{equation}
left{
egin{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}
ight .
end{equation}
]
对于（3）式中的第一式利用龙哥库塔方法，有：
[egin{equation}
left{
egin{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}
ight .
end{equation}
]
同理，对于(3)中的第二式：
[egin{equation}
left{
egin{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}
ight .
end{equation}
]
对于(3)中的第三式：
[egin{equation}
left{
egin{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}
ight .
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/
u 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