常微分方程数值解：欧拉公式

算法原理

对于常微分方程初值问题

在求解区间[a,b]上作等距分割的剖分，步长，记。用数值微商的方法，即用差商近似微商数值求解常微分方程。

用向前差商近似

做出y(x)的在x=x0处的一阶向前差商式： 

又，于是得到

而y(x1)的近似值y1可按

 或 

求得。类似地，由

 以及 

得到计算近似值的向前欧拉公式：

由差商（差分）得到的上述方程称为差分方程。

由yn直接算出yn+1值的计算格式称为显式格式，向前欧拉公式是显式格式。

算法流程

算法代码

[cpp]  view plain copy

1. 　　//欧拉公式代码  
2. #include<iostream>  
3. #include<string>  
4. #include<vector>  
5. using namespace std;  
6.   
7. double f(double x,double y){  
8.     return -50*y;  
9. }  
10.   
11. vector<double> Euler(double x0,double y0,double h,int N){  
12.     vector<double> Y(N,0);  
13.     double x=x0;  
14.     Y[0]=y0;  
15.     for(int n=1;n<N;n++){  
16.         Y[n] = Y[n-1] + h*f(x,Y[n-1]);  
17.         x += h;  
18.     }  
19.     return Y;  
20. }  
21.   
22. int main(){  
23.     char a='n';  
24.     do{  
25.         cout<<"请输入步长h和要计算的函数值的个数N： "<<endl;  
26.         double h;  
27.         int N;  
28.         cin>>h>>N;  
29.         cout<<"请输入要初始函数点（x0，y0）："<<endl;  
30.         double x0;  
31.         double y0;  
32.         cin>>x0>>y0;  
33.         vector<double> Y=Euler(x0,y0,h,N+1);  
34.         cout<<"欧拉格式计算结果为：  "<<endl;  
35.         for(int i=0;i<N+1;i++){  
36.             cout<<x0+i*h<<"     "<<Y[i]<<endl;  
37.         }  
38.         cout<<"是否要继续？(y/n)"<<endl;  
39.         cin>>a;  
40.     }while(a=='y');  
41.         return 0;  
42. }  
43.   
44.   
45.   
46. //改进的欧拉公式  
47. #include<iostream>  
48. #include<string>  
49. #include<vector>  
50. using namespace std;  
51.   
52. double f(double x,double y){  
53.     return y-(2*x)/y;  
54. }  
55.   
56. vector<double> ImprovedEuler(double x0,double y0,double h,int N){  
57.     vector<double> Y(N,0);  
58.     Y[0]=y0;  
59.     double x=x0;  
60.     double p=0;  
61.     double c=0;  
62.     for(int n=1;n<N;n++){  
63.         p=Y[n-1]+h*f(x,Y[n-1]);  
64.         x +=h;  
65.         c=Y[n-1]+h*f(x,p);  
66.         Y[n]=(p+c)/2;  
67.     }  
68.     return Y;  
69. }  
70.   
71. int main(){  
72.     char a='n';  
73.     do{  
74.         cout<<"请输入步长h和要计算的函数值的个数N： "<<endl;  
75.         double h;  
76.         int N;  
77.         cin>>h>>N;  
78.         cout<<"请输入要初始函数点（x0，y0）："<<endl;  
79.         double x0;  
80.         double y0;  
81.         cin>>x0>>y0;  
82.         vector<double> Y=ImprovedEuler(x0,y0,h,N+1);  
83.         cout<<"欧拉格式计算结果为：  "<<endl;  
84.         for(int i=0;i<N+1;i++){  
85.             cout<<x0+i*h<<"     "<<Y[i]<<endl;  
86.         }  
87.         cout<<"是否要继续？(y/n)"<<endl;  
88.         cin>>a;  
89.     }while(a=='y');  
90.         return 0;  
91. }  

实验过程原始记录

（1）分别取h=0.05，N=10；h=0.025，N=20；h=0.01，N=50，用显式欧拉方法求解微分方程初值 问题：y’=-50y，y(0)=10

h=0.05，N=10

h=0.025，N=20

h=0.01，N=50

（2）用改进的欧拉格式计算下列一阶常微分方程初值问题

其解析解为： 

实验结果及分析

1、欧拉公式用以求解常微分方程中的定解问题
2、可以看出，欧拉公式的精度很低，对于不同的步长求得相同点处的值差距可能很大；而且计算中的误差会累计。但显式欧拉公式取向前差商作为平均斜率，计算简单，且利于编写计算机程序，所以对于一些简单函数仍有很大的价值。
3、改进的欧拉公式是欧拉方法和梯形方法的综合，也是一种显式算法，计算简单，利于编写程序，与欧拉公式相比大大提高了精度。

（转载请注明作者和出处：http://blog.csdn.net/xiaowei_cqu 未经允许请勿用于商业用途）