欧拉方法
y n + 1 = y n + h f ( x n , y n ) , n = 0 , 1 , 2 , . . . y_{n+1}=y_{n}+hf(x_{n},y_{n}),n=0,1,2,... yn+1=yn+hf(xn,yn),n=0,1,2,...
是一种显式算法,计算量小,但精度很低。梯形方法
y n + 1 = y n + h 2 [ f ( x n , y n ) + f ( x n + 1 , y n + 1 ) ] y_{n+1}=y_{n}+\frac{h}{2}[f(x_{n},y_{n})+f(x_{n+1},y_{n+1})] yn+1=yn+2h[f(xn,yn)+f(xn+1,yn+1)]
虽然提高了精度,但它是一种隐式算法,需要借助于迭代过程求解,计算量大。
为此,我们将综合应用这两种方法,先用欧拉格式求得一个初步的近似值,记为 y ‾ n + 1 \overline{y}_{n+1} yn+1,称之为预报值;预报值 y ‾ n + 1 \overline{y}_{n+1} yn+1 的精度不高,我们用它替代梯形公式
y n + 1 = y n + h 2 [ f ( x n , y n ) + f ( x n + 1 , y n + 1 ) ] y_{n+1}=y_{n}+\frac{h}{2}[f(x_{n},y_{n})+f(x_{n+1},y_{n+1})] yn+1=yn+2h[f(xn,yn)+f(xn+1,yn+1)]
右端的 y ‾ n + 1 \overline{y}_{n+1} yn+1 再直接计算,得到校正值 y n + 1 y_{n+1} yn+1。这样建立的预报 - 校正系统
{ 预报: y ‾ n + 1 = y n + h f ( x n , y n ) 校正: y n + 1 = y n + h 2 [ f ( x n , y n ) + f ( x n + 1 , y ‾ n + 1 ) ] \begin{cases} 预报:\overline{y}_{n+1}=y_{n}+hf(x_{n},y_{n})\\ 校正:y_{n+1}=y_{n}+\frac{h}{2}[f(x_{n},y_{n})+f(x_{n+1},\overline{y}_{n+1})] \end{cases} {预报:yn+1=yn+hf(xn,yn)校正:yn+1=yn+2h[f(xn,yn)+f(xn+1,yn+1)]
称作改进的欧拉格式。这是一种一步显式格式,它可表为如下嵌套格式
y n + 1 = y n + h 2 [ f ( x n , y n ) + f ( x n + 1 , y n + h f ( x n , y n ) ) ] y_{n+1}=y_{n}+\frac{h}{2}[f(x_{n},y_{n})+f(x_{n+1},y_{n}+hf(x_{n},y_{n}))] yn+1=yn+2h[f(xn,yn)+f(xn+1,yn+hf(xn,yn))]
或表为下列平均化形式
{ y p = y n + h f ( x n , y n ) y c = y n + h f ( x n + 1 , y p ) y n + 1 = 1 2 ( y p + y c ) \begin{cases} y_{p}=y_{n}+hf(x_{n},y_{n})\\ y_{c}=y_{n}+hf(x_{n+1},y_{p})\\ y_{n+1}=\frac{1}{2}(y_{p}+y_{c}) \end{cases} ⎩ ⎨ ⎧yp=yn+hf(xn,yn)yc=yn+hf(xn+1,yp)yn+1=21(yp+yc)
图 3-2 描述了改进的欧拉方法,其中 h 为步长,N 为步数,x0,y0 为 “老值”,即前一步的近似解,x1,y1 为 “新值”,即该步计算结果。
例:用改进欧拉方法求解以下常微分方程的初值问题:
{ y ′ = y − 2 x y ( 0 < x < 1 ) y ( 0 ) = 1 \begin{cases} y'=y-\frac{2x}{y} (0 < x < 1)\\ y(0)=1 \end{cases} {y′=y−y2x(0<x<1)y(0)=1
解:求解初值问题的改进欧拉格式具有以下形式:
{ y p = y n + h ( y n − 2 x n y n ) y c = y n + h ( y p − 2 x n + 1 y p ) y n + 1 = 1 2 ( y p + y c ) \begin{cases} y_{p}=y_{n}+h(y_{n}-\frac{2x_{n}}{y_{n}})\\ y_{c}=y_{n}+h(y_{p}-\frac{2x_{n+1}}{y_{p}})\\ y_{n+1}=\frac{1}{2}(y_{p}+y_{c}) \end{cases} ⎩ ⎨ ⎧yp=yn+h(yn−yn2xn)yc=yn+h(yp−yp2xn+1)yn+1=21(yp+yc)
取 h = 0.1,计算结果如下表所示,同欧拉方法的计算结果比较,改进欧拉方法明显地改善了精度。
xn | yn | y(xn) |
---|---|---|
0.1 | 1.0959 | 1.0954 |
0.2 | 1.1841 | 1.1832 |
0.3 | 1.2662 | 1.2649 |
0.4 | 1.3434 | 1.3416 |
0.5 | 1.4164 | 1.4142 |
0.6 | 1.4860 | 1.4832 |
0.7 | 1.5525 | 1.5492 |
0.8 | 1.6165 | 1.6125 |
0.9 | 1.6782 | 1.6733 |
1.0 | 1.7379 | 1.7321 |
运行示例:
程序源码:
#include <iostream>
#include <iomanip>
#include <math.h>
using namespace std;
/**
* 欧拉格式中的 x,y 的函数关系式,即 f(xn,yn)
*/
double f(double x, double y)
{
return y - 2 * x / y;
}
/**
* 实际函数解析式
*/
double f1(double x)
{
return sqrt(1 + 2 * x);
}
int main(void)
{
double x0, y0;
cout << "请输入初值:";
cin >> x0 >> y0;
double h;
cout << "请输入步长:";
cin >> h;
int N;
cout << "请输入步数:";
cin >> N;
cout << "\t" << setw(10) << "xn"
<< "\t" << setw(10) << "yn"
<< "\t" << setw(10) << "y(xn)" << endl;
for (int i = 1; i <= N; i++)
{
// 离散点
double x1 = x0 + h;
// 根据欧拉格式求出的值,预报
double yp = y0 + h * f(x0, y0);
// 校正
double yc = y0 + h * f(x1, yp);
double y1 = (yp + yc) / 2;
// 输出本次步进后的离散点数据
cout << i << "\t" << setw(10) << x1 << "\t" << setw(10) << y1 << "\t" << setw(10) << f1(x1) << endl;
x0 = x1;
y0 = y1;
}
return 0;
}