微分方程-第一章习题5实验报告

第一章习题5实验报告一、实验目的通过用Euler 法、中点格式、预报－校正格式、Adams 预报—修正格式等算法求解一阶常微分方程初值问题000(,),()y f x y x x X y x y '=<<⎧⎪⎨=⎪⎩ 的数值解，经典四级四阶R-K 法是高精度单步法，四阶Adams 预估－校正法是线性多步法，这两个算法较实验一算法精度要高。

通过本次实验，掌握R-K 法的编程要领，掌握多步法初值的准备方法，深刻领会微分方程数值解的实质，体会单步法和线性多步法各自的优缺点，熟练掌握各算法的计算机实现过程，并能从理论及实验结果分析各种算法的优缺点。

二、实验内容用Euler 法、中点格式、预报－校正格式、经典四级四阶R-K 格式、Adams 预报—修正格式算法求解一阶常微分方程初值问题22d ,(1,2]22(1)1xy x x dy x y y ⎧=+∈⎪⎨⎪=⎩的数值解，其精确解为y = 三、算法1 Euler 法100(,)(),1,2,,j j j j y y hf x y y y x j N+=+⎧⎪⎨==⎪⎩ 其中0X x N h -⎡⎤=⎢⎥⎣⎦，h 为步长。

2 中点格式10(,(,))22(),1,2,,j j j j j j h h y y hf x y f x y y y x j N+⎧=+++⎪⎨⎪==⎩ 3 预报－校正格式1(0)(0)11100(,)((,)(,))2(),1,2,,j j j j j j j j j j y y hf x y h y y f x y f x y y y x j N++++⎧=+⎪⎪⎪=++⎨⎪⎪==⎪⎩4 四级四阶格式：11234121324300y (22)6(,)(,),0,1,2,,122(,)22(,)()j j j j j j j j j j h y K K K K K f x y h h K f x y K j N h h K f x y K K f x h y hK y x y +⎧=++++⎪⎪=⎪⎪⎪=++⎪=-⎨⎪=++⎪⎪⎪=++⎪=⎪⎩5 Adams 预报—修正格式 ()()()()⎪⎩⎪⎨⎧-+-++=--+-+=--+++---+修正格式预报格式2111132101519,924937595524n n n n n n n n n n n n n f f f y x f hy y f f f f h y y其中h 为步长。

四 程序一、Euler 法程序如下： #include <stdio.h>#define LEN 100int Euler(int a,int b,double h) {double y[LEN]={0}; double x[LEN]={0};int i=0;double num=0;num=((b-a)/h)+1;//节点个数//循环计算出各个节点的坐标 for(i=1;i<num;++i) { x[0]=a; x[i]=a+i*h;}//带公式计算各个节点的y的值for(i=0;i<num;i++){y[0]=1;y[i+1]=(y[i]/(2*x[i]))+(x[i]*x[i])/(2*y[i]);}//输出x和yprintf("\n ************************* Euler法求解方程******************************* \n");printf("\n x[n] y[n]\n\n");for(i=0;i<num;i++){printf(" %.4f %.4f\n",x[i],y[i]);}printf("\n ************************* Euler法求解方程******************************** \n\n");return 0;}}二、中点格式程序如下：int Mid(int a,int b,double h){double y[LEN]={0};double x[LEN]={0};int i=0;double num=0;num=((b-a)/h)+1;//节点个数//循环计算出各个节点的坐标for(i=0;i<num;++i){x[0]=a;x[i]=a+i*h;}//带公式计算各个节点的y的值for(i=0;i<num;i++){y[0]=1;y[i+1]=y[i]+h*f((x[i]+h/2),(y[i]+(h/2)*(f(x[i],y[i]))));}//输出x和yprintf("\n ************************* 中点格式求解方程******************************* \n");printf(" x[n] y[n]\n\n");for(i=0;i<num;i++){printf(" %.4f %.4f\n",x[i],y[i]);}printf("\n ************************* 中点格式求解方程******************************* \n\n");return 0;}三、预报——校正格式程序如下：int PreCorr(int a,int b,double h){double y1[LEN]={0};double y[LEN]={0};double x[LEN]={0};int i=0;double num=0;num=((b-a)/h)+1;//节点个数//循环计算出各个节点的坐标for(i=0;i<num;++i){x[0]=a;x[i]=a+i*h;}//带公式计算各个节点的y的值for(i=0;i<num;i++){y[0]=1;//初值为1y1[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],y1[i+1]));}//输出x和yprintf("\n *********************** 预报—校正格式求解方程*************************** \n");printf(" x[n] y[n]\n\n");for(i=0;i<num;i++){printf(" %.4f %.4f\n",x[i],y[i]);}printf("\n *********************** 预报—校正格式求解方程*************************** \n\n");return 0;}四、经典四级四阶R-K程序如下：int RungeKutta(int a,int b,double h){double y1[LEN]={0};double y[LEN]={0};double x[LEN]={0};int i=0;double num=0;num=((b-a)/h)+1;//节点个数//循环计算出各个节点的坐标for(i=0;i<num;++i){x[0]=a;x[i]=a+i*h;}//带公式计算各个节点的y的值for(i=0;i<num;i++){double K1=0;double K2=0;double K3=0;double K4=0;y[0]=1;//初值为1K1=f(x[i],y[i]);K2=f((x[i]+h/2),(y[i]+h*K1/2));K3=f((x[i]+h/2),(y[i]+h*K2/2));K4=f((x[i]+h),(y[i]+h*K3));y[i+1]=y[i]+(h/6)*(K1+2*K2+2*K3+K4);}//输出x和yprintf("\n ******************** 经典四级四阶R-K格式求解方程************************* \n");printf(" x[n] y[n]\n\n");for(i=0;i<num;i++){printf(" %.4f %.4f\n",x[i],y[i]);}printf("\n ******************** 经典四级四阶R-K格式求解方程************************* \n\n");return 0;}六、Adams程序如下：int Adams(int a,int b,double h){double y1[LEN] = {0};double y2[LEN] = {0};double y[LEN] = {0};double x[LEN] = {0};int i = 0;double num = 0;double tmp = 0;num = ((b-a)/h)+1;//节点个数//循环计算出各个节点的坐标for(i = 0; i < num; ++i){x[0] = a;x[i] = a + i*h;}for(i = 0; i < num; i++){y2[i] = sqrt(x[i]/2+(x[i]*x[i])/2);}y2[0] = 1;//带公式计算各个节点的y的值for(i = 0; i < 4; i++){double K1 = 0;double K2 = 0;double K3 = 0;double K4 = 0;y[0] = 1;//初值为1K1 = f(x[i],y[i]);K2 = f((x[i]+h/2),(y[i]+h*K1/2));K3 = f((x[i]+h/2),(y[i]+h*K2/2));K4 = f((x[i]+h),(y[i]+h*K3));y[i+1] = y[i]+(h/6)*(K1+2*K2+2*K3+K4);y1[i] = y[i];}for(i = 4;i < num;i++){y1[i+1] = y[i]+(h/24)*(55*f(x[i],y[i])-59*f(x[i-1],y[i-1])+37*f(x[i-2],y[i-2])-9*f(x[i-3],y[i-3]));y[i+1] = y[i]+(h/24)*(9*f(x[i+1],y1[i+1])+19*f(x[i],y[i])-5*f(x[i-1],y[i-1])+f(x[i-2],y[i-2]));}//输出x和yprintf("\n ***********************Adams预报—校正格式求解方程*************************** \n");printf(" x[n] y[n] y(精确解)\n\n");for(i = 0; i < num; i++){printf(" %.4f %.4f %.4f \n",x[i],y[i],y2[i]);}printf("\n **********************Adams预报—校正格式求解方程**************************** \n\n");return 0;}五、精确解程序如下：int ac_solve(int a,int b,double h){double y[LEN]={0};double x[LEN]={0};int i=0;double num=0;num=((b-a)/h)+1;//节点个数//循环计算出各个节点的坐标for(i=1;i<num;++i){x[0]=a;x[i]=a+i*h;}//带公式计算各个节点的y的值for(i=0;i<num;i++){y[i]=sqrt(x[i]/2+(x[i]*x[i])/2);}y[0]=1;//输出x和yprintf("\n **************************** 精确解***********************************\n");printf("\n x[n] y[n]\n\n");for(i=0;i<num;i++){printf(" %.4f %.4f\n",x[i],y[i]);}printf("\n **************************** 精确解***********************************\n\n");return 0;}调用前四个函数#include<stdio.h>#include<math.h>#include "Euler.h"#define LEN 100int main(){int a=1;//区间左端点int b=2;//区间右端点int case_flag=0;double h=0.1;//步长item();scanf("%d",&case_flag);printf("\n");switch(case_flag){case 0:{break;}case 1:{Euler(a,b,h);//Euler法break;}case 2:{Mid(a,b,h);//中点格式break;}case 3:{PreCorr(a,b,h);//预报—校正格式break;}case 4:{RungeKutta(a,b,h);//经典四级四阶R-K格式break;}case 5:{ac_solve(a,b,h);//精确解break;}default:{printf("没有您所选择的序号!!!\n");break;}}return 0;}调用Adams预报—校正格式主函数int main(void){int a = 1;//区间左端点int b = 2;//区间右端点double h = 0.1;//步长Adams(a,b,h);//Adams法return 0;}五、数值结果及分析运算结果可用下面表格表示：由以上表格可知：从表5.1和表5.2可以看出在四种方法中Euler格式计算的方程的解与精确解相差最大，即误差较大；中点格式，预报-校正格式和经典四级四阶R-K格式的误差较小，其中经典四级四阶R-K格式离精确解最近，所以，经典四级四阶R-K格式方法在前四种方法中求解方程的解的误差最小，效果最好。