求简单的常微分方程

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

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

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

　　参考：常微分方程组的数值算法解法 P_1~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})] $$

　　由于利用该公式进行求解时需要求解含有y_i+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(x_n+1)-y_n+1有多大？这里假定前一步得的结果y_n=y(x_n)是准确的。

　　　　(1)欧拉法的截断误差

　　　　　　将y(x_n+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(h²)，当h→0时，它与h²是同阶无穷小量。

　　　　(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} $$

　　　　　　将k₁、k₂代入预报校正公式得：

$$ \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} $$

　　略去截断误差，并用近似值y_n^(k)代替真值y^(k)(x_n)则得

$$ 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} $$

龙格-库塔法：

　　

转载于:https://www.cnblogs.com/darkchii/p/9175716.html