求简单的常微分方程

Posted 2020-11-17 darkchii

求常微分方程的原理(懒得重新打一遍。。于是把我知乎上的一个相关回答搬过来)：

这里介绍一种方法，叫欧拉法，比如，形如：

$$ left{ egin{gathered} quad frac{dy}{dx}=f(x,y) \ y(x_0)=y_0 end{gathered} ight. $$

的一阶微分方程（注：用数值方法求解时，默认有解且唯一）。

通过初始条件： $$ y(x_0)=y_0 $$

可以计算出： $$ y‘(x_0)=f(x_0,y_0) $$

假设  充分小，则近似的有：

$$ 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为步长 $$

参考：常微分方程组的数值算法解法 P_1~2

 

例题：求一阶微分方程：

$$ left{ egin{gathered} frac{dy}{dx}=y \ y(0)=1 end{gathered} ight. \ 步长h=0.01，x=0.10时的值 $$

下面是代码实现：

#include<bits/stdc++.h>
int main()
{
    // 欧拉法：yi+1 = yi + h * f(x,y)
    // 
    double h = 0.001, x = 0.10, y = 1.00;
    for(double i = 0.00; i <= x; i += h)
        y += h*y;
    printf("%.4f
", 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.

看起来还是不错的，但精度还是不够，一些特殊方程，或许需要更高的精度，则会导致计算效率非常差。。。

使用龙格-库塔法提高精度：

ps:明天再填。。。