求常微分方程的原理(懒得重新打一遍。。于是把我知乎上的一个相关回答搬过来)：
这里先介绍一种方法，叫欧拉法，比如，形如：
$$ left{ egin{aligned} frac{dy}{dx} &= f(x,y) \ y(x_0) &= y_0 end{aligned} ight. $$
的一阶微分方程（注：用数值方法求解时，默认有解且唯一）。
通过初始条件： $$ 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{ egin{aligned} frac{dy}{dx} &= y \ y(0) &= 1 end{aligned} ight. \ 步长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
", 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 { egin{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} ight. $$
通常将上式称为预报校正公式，其中第一式叫预报公式，第二式叫校正公式。这个公式还可写为：
$$ left { egin{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} ight. $$
利用这个方法测试一下例题1：
Output:
[darkchii@localhost 文档]$ ./test Euler:1.1157 ReEuler:1.1163
emmm...貌似我的写法有问题？
例题2：求一阶微分方程
$$ left{ egin{aligned} frac{dy}{dx} &= y-frac{2x}{y} \ y(0) &= 1 end{aligned} ight. $$
$$ 当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
", 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)用泰勒展开：
$$ egin{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)改进欧拉法的截断误差
以迭代一次的预报校正公式来说明，因为
$$ egin{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) \ &= hfegin{bmatrix} x_n + h, y(x_n) + k_1 end{bmatrix} \ &= h egin{Bmatrix} fegin{bmatrix}x_n, y(x_n)end{bmatrix} + h frac{partial}{partial x} fegin{bmatrix} x_n, y(x_n)end{bmatrix} + k_1 frac{partial}{partial y} fegin{bmatrix} x_n, y(x_n)end{bmatrix} + cdot cdot cdot end{Bmatrix} \ &=hfegin{bmatrix} x_n, y(x_n)end{bmatrix} + h^2 egin{Bmatrix} frac{partial}{partial x} fegin{bmatrix} x_n, y(x_n)end{bmatrix} + y'(x_n) frac{partial}{partial y} fegin{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代入预报校正公式得：
$$ egin{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{ egin{aligned} frac{dy}{dx} &= f(x,y) \ y(x_0) &= y_0 end{aligned} ight. $$
的精确解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 $$
其中截断误差为
$$ egin{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} $$
龙格-库塔法：