求常微分方程的原理(懒得重新打一遍。。于是把我知乎上的一个相关回答搬过来):
这里先介绍一种方法,叫欧拉法,比如,形如:
$$ \left\{ \begin{aligned} \frac{dy}{dx} &= f(x,y) \\ y(x_0) &= y_0 \end{aligned} \right. $$
的一阶微分方程(注:用数值方法求解时,默认有解且唯一)。
通过初始条件: $$ y(x_0)=y_0 $$
可以计算出: $$ y'(x_0)=f(x_0,y_0) $$
假设 $ x_1 - x_0 = h $ 充分小,则近似的有:
$$ \frac{y(x_1)-y(x_0)}{h} \approx y'(x_0)=f(x_0,y_0) \quad $$
记 $$ y_i=y(x_i) \quad i=0,1,...,n $$
取 $$ y_1=y_0+hf(x_0,y_0) $$
同样的方法,计算出
$$ y_2=y_1+hf(x_1,y_1) $$
于是得到递推公式:
$$ y_{i+1}=y_i+hf(x_i,y_i) ,h为步长 $$
参考:常微分方程组的数值算法解法 P1~2
p.s.可看到实际上与牛顿切线法原理类似。。。牛顿法:
$$ x_{i+1}=x_i-\frac{f(x_i)}{f'(x_i)} $$
例题1:求一阶微分方程
$$ \left\{ \begin{aligned} \frac{dy}{dx} &= y \\ y(0) &= 1 \end{aligned} \right. \\ 步长h=0.01,x=0.10时的值 $$
下面是代码实现:
#include<bits/stdc++.h>
int main()
{
// 欧拉法:yi+1 = yi + h * f(x,y)
//
double h = 0.01, x = 0.10, y = 1.00;
for(double i = 0.00; i < x; i += h)
y += h*y;
printf("%.4f\n", y);
return 0;
}
Output:
1.1157
Process returned 0 (0x0) execution time : 0.420 s
Press any key to continue.
精确结果:
e^0.1
1.1051709180756
当h=0.001时:
1.1051
Process returned 0 (0x0) execution time : 0.401 s
Press any key to continue.
看起来还是不错的,但精度还是不够,一些方程还需要更高的精度且x越大则会导致计算效率非常差。。。
改进欧拉法:
考虑下面积分方程:
$$ y(x)=y_0+ \int _{x_0}^x f(t,y(t))dt $$
等式右边第二项利用积分中值定理有:
$$ \int _{x_0}^x f(t,y(t))dt \approx f(x_0,y(x_0))(x-x_0) $$
$$ 当 x=x_1时,有\int _{x_0}^{x_1} f(t,y(t))dt \approx f(x_0,y(x_0))(x_1-x_0) $$
$$ 取x_1-x_0=h,于是有y1=y0+hf(x_0,y(x_0)) $$
这个形式与前面欧拉法递推公式相同,所以我们可以理解为什么欧拉法精度上会有不足,因为等式右边相当于利用黎曼求和法进行了函数的矩形面积求和,所以我们可以考虑利用梯形法求和提高精度,利用梯形公式有:
$$ \int _{x_0}^{x_1} f(t,y(t))dt \approx \frac{1}{2}[f(x_0,y(x_0))+f(x_1,y(x_1))](x_1-x_0) $$
$$ 取x_1-x_0=h,有: $$
$$ y_1=y_0+\frac{h}{2}[(f(x_0,y_0)+f(x_1,y_1)] $$
于是递推公式变为:
$$ y_{i+1}=y_i+\frac{h}{2}[(f(x_i,y_i)+f(x_{i+1},y_{i+1})] $$
由于利用该公式进行求解时需要求解含有yi+1的方程,这常常是困难的,所以实际计算时,常常与欧拉法结合使用:
$$ \left \{ \begin{aligned} y_{i+1}^{(0)} &= y_i+hf(x_i,y_i) \\ \quad y_{i+1}^{(k+1)} &= y_i+\frac{h}{2}[(f(x_i,y_i)+f(x_{i+1},y_{i+1}^{(k)})] \quad k=0,1,2... \end{aligned} \right. $$
通常将上式称为预报校正公式,其中第一式叫预报公式,第二式叫校正公式。这个公式还可写为:
$$ \left \{ \begin{aligned} y(x_{n+1}) &= y_n + \frac{1}{2}k_1 + \frac{1}{2}k_2 \\ k_1 &= hf(x_n, y_n) \\ k2 &= hf(x_n + h, y_n + k_1) \end{aligned} \right. $$
利用这个方法测试一下例题1:
Output:
[darkchii@localhost 文档]$ ./test
Euler:1.1157 ReEuler:1.1163
emmm...貌似我的写法有问题?
例题2:求一阶微分方程
$$ \left\{ \begin{aligned} \frac{dy}{dx} &= y-\frac{2x}{y} \\ y(0) &= 1 \end{aligned} \right. $$
$$ 当h=0.1,x=1.5时的值 $$
代码实现:
#include<bits/stdc++.h>
int main()
{
double h = 0.1, x = 1.50, y = 1.00, result = y;
// y_{i+1} = y_i + h * f(x_i,y_i)
// y_{i+1} = y_i + h/2 * (f(x_i,y_i) + f(x_{i+1},y_{i+1}))
for(double i = 0.00; i < x; i += h)
{
y += h*(y-((2*i)/y));
result += (h/2.0)*(result-((2*i)/result) + y-((2*(i+h))/y));
}
printf("Euler:%.4f ReEuler:%.4f\n", y, result);
return 0;
}
Output:
[darkchii@localhost 文档]$ ./test
Euler:2.1201 ReEuler:2.0819
改进的欧拉法结果与书上例题给出的结果没对上,我的写法可能有问题。。。但按书上的意思是改进欧拉法只使用第一次欧拉法迭代一次的结果,但我觉得只用一次貌似没有意义...
现在来考察两个公式的截断误差:即y(xn+1)-yn+1有多大?这里假定前一步得的结果yn=y(xn)是准确的。
(1)欧拉法的截断误差
将y(xn+1)用泰勒展开:
$$ \begin{aligned} y(x_{n+1}) &= y(x_n + h) \\ &= y(x_n) + hy'(x_n) + \frac{h^2}{2!}y''(x_n) + \cdot \cdot \cdot \quad (1) \end{aligned} $$
由欧拉法得:
$$ y_{n+1} = y_n + hf(x_n, y_n) = y(x_n) + hy'(x_n) \quad (2) $$
将(1)、(2)两式相减有:
$$ y(x_{n+1}) - y_{n+1} = \frac{h^2}{2} y''(x_n) + \cdot \cdot \cdot = O(h^2) $$
即欧拉法的截断误差为O(h2),当h→0时,它与h2是同阶无穷小量。
(2)改进欧拉法的截断误差
以迭代一次的预报校正公式来说明,因为
$$ \begin{aligned} k_1 &= hf(x_n, y_n) = hf(x_n, y(x_n)) \\ &= hy'(x_n) \\ k_2 &= hf(x_n + h, y_n + k_1) \\ &= hf\begin{bmatrix} x_n + h, y(x_n) + k_1 \end{bmatrix} \\ &= h \begin{Bmatrix} f\begin{bmatrix}x_n, y(x_n)\end{bmatrix} + h \frac{\partial}{\partial x} f\begin{bmatrix} x_n, y(x_n)\end{bmatrix} + k_1 \frac{\partial}{\partial y} f\begin{bmatrix} x_n, y(x_n)\end{bmatrix} + \cdot \cdot \cdot \end{Bmatrix} \\ &=hf\begin{bmatrix} x_n, y(x_n)\end{bmatrix} + h^2 \begin{Bmatrix} \frac{\partial}{\partial x} f\begin{bmatrix} x_n, y(x_n)\end{bmatrix} + y'(x_n) \frac{\partial}{\partial y} f\begin{bmatrix}x_n, y(x_n)\end{bmatrix} + \cdot \cdot \cdot \end{Bmatrix} \\ &= hy'(x_n) + h^2 y''(x_n) + \cdot \cdot \cdot \end{aligned} $$
将k1、k2代入预报校正公式得:
$$ \begin{aligned} y_{n+1} &= y_n + \frac{h}{2} y'(x_n) + \frac{h}{2} y'(x_n) + \frac{h^2}{2} y''(x_n) + \cdot \cdot \cdot \\ &= y(x_n) + hy'(x_n) + \frac{h^2}{2}y''(x_n) + \cdot \cdot \cdot \quad (3) \end{aligned} $$
再将(1)、(3)两式相减有:
$$ y(x_{n+1}) - y_{n+1} = O(h^3) $$
即改进欧拉法比欧拉法的阶提高了。
泰勒级数法:
如果初值问题
$$ \left\{ \begin{aligned} \frac{dy}{dx} &= f(x,y) \\ y(x_0) &= y_0 \end{aligned} \right. $$
的精确解y(x)在[x,x]上存在k+1阶导数且连续,那么由泰勒公式
$$ y(x_{n+1}) = y(x_n) + hy'(x_n) + \cdot \cdot \cdot + \frac{h^k}{k!} y^{k}(x_n) + R_k $$
其中截断误差为
$$ \begin{aligned} R_k= \frac{h^{k+1}}{(k+1)!} y^{k+1}(\xi) = O(h^{k+1}) \\ x_n < \xi < x_{n+1} \end{aligned} $$
略去截断误差,并用近似值yn(k)代替真值y(k)(xn)则得
$$ y_{n+1} = y_n + hy'_n + \frac{h^2}{2}y''_n + \cdot \cdot \cdot + \frac{h^k}{k!}y_n ^{(k)} $$
做数值计算时一般取k=3阶就足够了,这时的截断误差为
$$ R_3= \frac{h^{4}}{(4)!} y^{4}(\xi) = O(h^{4}) \\ x_n < \xi < x_{n+1} $$
龙格-库塔法: