今天学习常微分方程的数值解，学习视频指路：微分方程数值解，感谢本视频让我突击学会了数值解（教材真的看不懂www）。本文作为该视频的学习笔记，加入了个人的一些理解，同时参考了王高雄老师的《常微分方程》第四版，如有理解不当地方还请各位指正。
文章目录
一、 欧拉方法（Euler法）
01 Euler法基本格式
对于微分方程的初值问题 { d y d x = f ( x , y ) y ( x 0 ) = y 0 \begin{cases}\frac{dy}{dx}=f(x,y)\\ y(x_0)=y_0\end{cases} {dxdy=f(x,y)y(x0)=y0数值方法本质上就是希望能够找出一系列点去逼近 y = f ( x ) y=f(x) y=f(x)这条曲线。
假设这个微分方程的解曲线在
[
0
,
T
]
[0,T]
[0,T]上有定义，我们将区间划分为
n
n
n个区间，即
n
+
1
n+1
n+1个节点，并令步长
h
=
T
n
h=\frac{T}{n}
h=nT,记节点为
t
k
t_k
tk,满足
0
<
t
0
<
t
1
<
.
.
.
<
t
n
=
T
0<t_0<t_1<...<t_n=T
0<t0<t1<...<tn=T
为了方便推导符号的统一性（其实用
y
y
y也是可以的），我们记初值问题为
{
u
(
t
0
)
=
u
0
u
′
(
t
0
)
=
f
(
t
0
,
u
0
)
\begin{cases}u(t_{0})=u_{0}\\u'(t_{0})=f(t_0,u_{0})\end{cases}
{u(t0)=u0u′(t0)=f(t0,u0)在下一个节点处进行泰勒展开有
u
(
t
1
)
=
u
(
t
0
+
h
)
=
u
(
t
0
)
+
u
′
(
t
0
)
h
+
O
(
h
2
)
u(t_1)=u(t_0+h)=u(t_0)+u'(t_0)h+O(h^2)
u(t1)=u(t0+h)=u(t0)+u′(t0)h+O(h2)忽略二阶小量，令一个新的符号表示近似值
u
1
=
u
(
t
0
)
+
u
′
(
t
0
)
h
u_1=u(t_0)+u'(t_0)h
u1=u(t0)+u′(t0)h由于在递推过程中，精确解本身就是难以得到的，所有的逼近过程都是用近似值
u
1
=
u
0
+
f
(
t
0
,
u
0
)
h
≈
u
(
t
1
)
u_1=u_0+f(t_0,u_0)h\approx u(t_1)
u1=u0+f(t0,u0)h≈u(t1)以此类推，得到递推公式
u
k
+
1
=
u
k
+
f
(
t
k
,
u
k
)
h
(
k
=
0
,
1
,
.
.
.
,
n
−
1
)
(1-1)
u_{k+1}=u_k+f(t_{k},u_{k})h\quad (k=0,1,...,n-1)\tag{1-1}
uk+1=uk+f(tk,uk)h(k=0,1,...,n−1)(1-1)其几何意义为通过递推式子，我们可以得到各个节点的函数值，连接可以得到一个折线（其实也可以说是用差分近似微分），每一点都是取其上一点的切线近似而成（切线也不是精确的，是近似的切线），就是以直代曲的思想，这个折线就是解曲线的近似曲线，如果区间划分的更细致（
n
n
n更大）时，这个折线会逐渐趋近于解曲线，达到更好的近似效果。
02 差分形式拓展
如果我们用差分代替微分的方式来理解公式(1-1)的话，可以将其改写为
u
′
(
t
k
)
≈
u
k
+
1
−
u
k
h
=
f
(
t
k
,
u
k
)
u'(t_k)\approx\frac{u_{k+1}-u_k}{h}=f(t_k,u_k)
u′(tk)≈huk+1−uk=f(tk,uk)这是向后差分，也称为显式欧拉法。当然肯定会有同学有问号，为什么不能用
u
k
−
u
k
−
1
u_{k}-u_{k-1}
uk−uk−1来近似，事实上是可以的，如下
u
′
(
t
k
)
≈
u
k
−
u
k
−
1
h
=
f
(
t
k
,
u
k
)
u'(t_k)\approx\frac{u_{k}-u_{k-1}}{h}=f(t_k,u_k)
u′(tk)≈huk−uk−1=f(tk,uk)这是向前差分，
u
k
u_k
uk难以直接求出，该方法称为隐式欧拉法。
03 Euler法的改进
在数值方法中经常使用积分近似：
u
(
t
)
=
u
(
t
0
)
+
∫
t
0
t
f
(
τ
,
u
(
τ
)
)
d
τ
u(t)=u(t_0)+\int_{t_0}^tf(\tau,u(\tau))d\tau
u(t)=u(t0)+∫t0tf(τ,u(τ))dτ由于积分是要选定积分区间的，当采用近似方法时，在一个非常小的划分区间内，我们将被积分的导数看作是常数，被积分的导数可以选择左端点的导数，也可以选择右端点的导数，由此衍生出三种近似方法。
（1）左矩形公式
以左边端点的导数进行近似
u
1
=
u
0
+
f
(
t
0
,
u
0
)
h
u_1=u_0+f(t_0,u_0)h
u1=u0+f(t0,u0)h该方法实际上就是上面所述的显式欧拉法。
（2）右矩形公式
以右边端点的导数进行近似
u
1
=
u
0
+
f
(
t
1
,
u
1
)
h
u_1=u_0+f(t_1,u_1)h
u1=u0+f(t1,u1)h该方法实际上就是上面所述的隐式欧拉法。
（3）梯形公式
为了更好的近似效果，我们希望将两个端点都考虑进去。以两个端点的导数算术平均数作为切线斜率，得到梯形公式
u
1
=
u
0
+
f
(
t
0
,
u
0
)
+
f
(
t
1
,
u
1
)
2
h
u_1=u_0+\frac{f(t_0,u_0)+f(t_1,u_1)}{2}h
u1=u0+2f(t0,u0)+f(t1,u1)h利用梯形公式就可以得到该情况下欧拉迭代公式
u
k
+
1
=
u
k
+
f
(
t
k
,
u
k
)
+
f
(
t
k
+
1
,
u
k
+
1
)
2
h
u_{k+1}=u_k+\frac{f(t_k,u_k)+f(t_{k+1},u_{k+1})}{2}h
uk+1=uk+2f(tk,uk)+f(tk+1,uk+1)h这实际上也是一种隐式的欧拉法，因为
u
k
+
1
u_{k+1}
uk+1在等式的左右两端都出现了，因此
u
k
+
1
u_{k+1}
uk+1实际是不能直接求出来的。为了解决这个问题，我们可以先用常规的欧拉法进行预测，得到一个近似的
u
k
+
1
u_{k+1}
uk+1，再用梯形公式进行修正，得到改进的欧拉法。
u
ˉ
k
+
1
=
u
k
+
f
(
t
k
,
u
k
)
h
u
k
+
1
=
u
k
+
f
(
t
k
,
u
k
)
+
f
(
t
k
+
1
,
u
ˉ
k
+
1
)
2
h
{\bar u_{k+1}}=u_k+f(t_k,u_k)h\\ u_{k+1}=u_k+\frac{f(t_k,u_k)+f(t_{k+1},\bar{u}_{k+1})}{2}h
uˉk+1=uk+f(tk,uk)huk+1=uk+2f(tk,uk)+f(tk+1,uˉk+1)h也就是我们将
u
ˉ
k
+
1
\bar{u}_{k+1}
uˉk+1作为一个中间量，代替了
f
(
t
k
+
1
,
u
k
+
1
)
f(t_{k+1},u_{k+1})
f(tk+1,uk+1)中的
u
k
+
1
u_{k+1}
uk+1，避免了通过隐式格式求解
u
k
+
1
u_{k+1}
uk+1的麻烦
04 实际运行效果
下图是在原视频里面截取的，运行软件为MATLAB
上面的图片分别是步长为0.05和0.01的近似图，可以看出当区间划分的越细时，近似效果越好，越接近真实曲线。
二、 龙格—库塔方法
三、 数值方法的评价
我们定义局部截断误差为
u
(
t
k
)
−
u
k
=
O
(
h
2
)
u(t_k)-u_{k}=O(h^2)
u(tk)−uk=O(h2)对于一个数值方法，我们最关心的一定是其误差如何，如果截断误差越小，意味着该方法的效果越好，我们可以看到阶段误差一般是某个量的高阶无穷小，对此我们定义，若某种数值方法的局部截断误差为步长
h
h
h的
O
(
h
p
+
1
)
O(h^{p+1})
O(hp+1)时，称该方法有
p
p
p阶精度。
欧拉法推导中使用了泰勒展开，从其展式中可以看出欧拉法有1阶精度。
扫一扫
1385
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
