数值分析常微分数值解的求法C语言

本科生课程设计报告实习课程数值分析
学院名称管理科学学院
专业名称信息与计算科学
学生姓名
学生学号
指导教师
实验地点
实验成绩
二〇一六年六月二〇一六年六
摘要
常微分方程数值解法是计算数学的一个分支.是解常微分方程各类定解问题的数值方法.现有的解析方法只能用于求解一些特殊类型的定解问题，实用上许多很有价值的常微分方程的解不能用初等函数来表示，常常需要求其数值解.所谓数值解，是指在求解区间内一系列离散点处给出真解的近似值.这就促成了数值方法的产生与发展.?
关键词：数值解法；常微分
1. 实验内容和要求
常微分方程初值问题 有精确解2()cos(2)x y x x e x -=+。

要求：分别取步长h = 0.1，0.01，0.001，采用改进的Euler 方法、4阶经典龙格-库塔R －K 方法和4阶Adams 预测-校正方法计算初值问题。

以表格形式列出10个等距节点上的计算值和精确值，并比较他们的计算精确度。

其中多步法需要的初值由经典R-K 法提供。

运算时，取足以表示计算精度的有效位。

2. 算法说明
2.1函数及变量说明
表1 函数及变量定义
1、欧拉方法：
1()()(,())i i i i y x y x hf x y x +=+ （i=0，1，2，3，......n-1）
(0)y a
= (其中a 为初值)
2、改进欧拉方法：
~
1~111()()(,())
()()[(,())(,())]
2
(0)i i i i i i i i i i y x y x hf x y x h
y x y x f x y x f x y x y a ++++=+=++=（i=0，1，2......n-1） (其中a 为初值)
3、经典K-R 方法： 1
1213
243()6
(,)(,)22(,)22
(,)
i i i i i i i i i i h y y K f x y h h
K f x y K h h K f x y K K f x h y hK +⎧
=+⎪⎪
=⎪⎪⎪=++⎨⎪
⎪
=++⎪⎪=++⎪⎩
4、4阶adams 预测-校正方法 预测： 校正：
Adsms 内插外插公式联合使用称为Adams 预测-校正系统，利用外插公式计算预测，用内插公式进行校正。

计算时需要注意以下两点： 1、外插公式为显式，内插公式为隐式。

故用内插外插公式时需要进行迭代。

2、这种预测-校正法是四步法，计算Yn+1时，不但用到前一步信息，而且要用到更前三步信息1-n f ，2-n f 3-n f ，因此它不是自动开始的，实际计算时必须借助某种单步法，用Runge-Kutta 格式为它提供初始值
3
21,,y y y ，依据局部截断误差公式得：
进一步将Adams 预测-校正系统加工成下列方案：
运用上述方案计算1n y +时，要先一步的信息n n n
c p y -,,y '
n 和更前一步的信息1-n y
因此在计算机之前必须给出初值1y 和11c -p ，1y 可用其他单步法来计算，11c -p 则一般令他为零。

3. 源程序
#include<stdio.h>
#include<stdlib.h>
#include<windows.h>
#include<math.h>
//微分方程
double f(double x,double y)
{
return (-y+cos(2*x)-2*sin(2*x)+2*x*exp(-x));
//return x*pow(y,-2)*2.0/3.0;
}
//原函数
double p_f(double x)
{
return x*x*exp(-x)+cos(2*x);
//return pow(1.0+pow(x,2),1.0/3.0);
}
//求精确解
double* Exact(double a,double b,double h)
{
int i;
double c = (b-a)/h;
double *y = new double [c+1];
for(i=0;i<=c;i++)
y[i]= p_f(a+i*h);
return y;
}
//欧拉方法
double* Euler(double a,double b,double h,double y0) {
int i;
double c = (b-a)/h;
double *y = new double [c+1];
y[0]=y0;
for(i=1;i<=c;i++)
{
y[i]= y[i-1]+ h* f(a+(i-1)*h,y[i-1]);
//printf("%f\n",f(a+(i-1)*h,y[i-1]));
}
return y;
}
//改进欧拉方法
double* Euler_Pro(double a,double b,double h,double y0) {
int i;
double c = (b-a)/h ,yb;
double *y = new double [c+1];
y[0]=y0;
for(i=1;i<=c;i++)
{
yb=y[i-1] +h* f(a+(i-1)*h,y[i-1]);
y[i] = y[i-1] + 0.5*h*( f(a+(i-1)*h,y[i-1]) + f(a+i*h,yb) );
}
return y;
}
//经典K-R方法
double* K_R(double a,double b,double h,double y0) {
double K1,K2,K3,K4,x;
int i;
double c = (b-a)/h ;
double *y = new double [c+1];
y[0]=y0;
for(i=1;i<=c;i++)
{
x=a+(i-1)*h;
K1=f(x,y[i-1]);
K2=f(x+0.5*h,y[i-1]+0.5*h*K1);
K3=f(x+0.5*h,y[i-1]+0.5*h*K2);
K4=f(x+h,y[i-1]+h*K3);
y[i]=y[i-1] + h*( K1 + 2*K2 + 2*K3 + K4)/6;
}
return y;
}
//4阶Adams预测-校正方法
double **Adams(double a,double b,double h,double y0)
{
int i;
double *y;
double count = (b-a)/h;
double dy[100]={0},x[4]={0};
double p_0=0,p_1=0,c=0;
double **r = new double *[count+1];//动态初始化，储存预测值和校值
for(i=0;i<=count;i++)
r[i] = new double [2];
if(r==NULL)
{
printf("内存分配失败\n");
exit(0);
}
for(i=0;i<count;i++)
{
r[i][0]=0;
r[i][1]=0;
}
y=K_R(a,b,h,y0);
for(i=0;i<4;i++)
{
x[i]=a+h*i;
dy[i]=f(x[i],y[i]);
r[i][0]=y[i];
}
i=3;
while(x[3]<b)
{
x[3] = x[3]+h;
p_1=y[3]+h*(55*dy[3]-59*dy[2]+37*dy[1]-9*dy[0])/24; //预估
c=p_1+251*(c-p_0)/270; //修正
r[i][0]=c;
c=f(x[3],c); //求f
c=y[3]+h*(9*c+19*dy[3]-5*dy[2]+dy[1])/24; //校正
y[3]=c-19*(c-p_1)/270; //修正
r[i][1] = y[3];
dy[0]=dy[1];
dy[1]=dy[2];
dy[2]=dy[3];
dy[3]=f(x[3],y[3]);
p_0=p_1;
i++;
}
return r;
}
void main()
{
int i , selet;
double a=0, b=2, h=0.001, y0=1;
while(1)
{
printf("请输入下限，上限，步长，初值：(用空格隔开)\n");
scanf("%lf %lf %lf %lf",&a,&b,&h,&y0);
double c = (b-a)/h;
double *yx,*y;
yx=Exact(a,b,h);
while(1)
{
printf("1、欧拉 2、改进欧拉 3、经典K-R4、4阶Adams 预测-校正 5、修改数据 6、退出\n");
scanf("%d",&selet);
system("cls");
if(selet==1)
{
y=Euler(a,b,h,y0);
printf("欧拉方法：\n");
}
if(selet==2)
{
y=Euler_Pro(a,b,h,y0);
printf("改进欧拉方法：\n");
}
if(selet==3)
{
y=K_R(a,b,h,y0);
printf("经典K-R方法(4阶龙格-库塔方法)：\n");
}
if(selet==4)
{
double **r;
r=Adams(a,b,h,y0);
printf("4阶Adams预测-校正方法:\n");
printf("x值预估值校正值误差\n");
printf("%0.3f %lf -
-\n",0,r[0][0]);
printf("%0.3f %lf -
-\n",0.1,r[1][0]);
printf("%0.3f %lf -
-\n",0.2,r[2][0]);
for(i=3;i<=10;i++)
printf("%0.3f %f %f %0.20f\n",a+i*h,r[i][0],r[i] [1],fabs(r[i][0]-r[i][1]));
}
if(selet==1||selet==2||selet==3)
{
printf("x值计算值准确值误差\n");
for(i=0;i<=10;i++)
printf("%0.3f %lf %lf %0.20f\n",(a+i)*h,y[i],yx[i] ,fabs(yx[i]-y[i]));
}
if(selet==5)
break;
if(selet==6)
exit(1);
continue;
}
}
}
4.实验结果
1、改进欧拉方法
2、经典K-R方法
3、4阶adams预测-校正法
5.对比分析
通过表格可以看出：
1，在步长相同的时候，结果准确度经典K-R>改进欧拉>欧拉方法>4阶adams预测校正方法
2，同一个方法的情况下，步长越小结果准确度越高；
参考文献：
[1] 朱建新，李有法，数值计算方法(第三版)，高等教育出版社；
[2]Adam预测-校正系统。