數(shù)值分析常微分?jǐn)?shù)值解的求法C語言

1、本科生課程設(shè)計報告實習(xí)課程數(shù)值分析學(xué)院名稱管理科學(xué)學(xué)院專業(yè)名稱信息與計算科學(xué)學(xué)生姓名學(xué)生學(xué)號指導(dǎo)教師實驗地點實驗成績二 一 六 年 六 月 二 一 六 年 六摘要常微分方程數(shù)值解法是計算數(shù)學(xué)的一個分支.是解常微分方程各類定解問題的數(shù)值方法.現(xiàn)有的解析方法只能用于求解一些特殊類型的定解問題，實用上許多很有價值的常微分方程的解不能用初等函數(shù)來表示，常常需要求其數(shù)值解.所謂數(shù)值解，是指在求解區(qū)間內(nèi)一系列離散點處給出真解的近似值.這就促成了數(shù)值方法的產(chǎn)生與發(fā)展.關(guān)鍵詞：數(shù)值解法；常微分1. 實驗內(nèi)容和要求常微分方程初值問題有精確解。要求：分別取步長h = 0.1，0.01，0.001，采用改進的Eul

2、er方法、4階經(jīng)典龍格-庫塔RK方法和4階Adams預(yù)測-校正方法計算初值問題。以表格形式列出10個等距節(jié)點上的計算值和精確值，并比較他們的計算精確度。其中多步法需要的初值由經(jīng)典R-K法提供。運算時，取足以表示計算精度的有效位。2. 算法說明2.1函數(shù)及變量說明表1 函數(shù)及變量定義函數(shù)/變量名類 型說 明fdouble要求的微分方程p_fdouble微分方程原函數(shù)Euler，Euler_ProK_R,Adams double歐拉，改進歐拉，經(jīng)典K_R，4階Adams預(yù)測-校正方法MainVoid主函數(shù)1、 歐拉方法： （i=0，1，2，3，.n-1） (其中a為初值)2、 改進歐拉方法：（i=

3、0，1，2.n-1） (其中a為初值)3、 經(jīng)典K-R方法： 4、4階adams預(yù)測-校正方法預(yù)測：校正： Adsms內(nèi)插外插公式聯(lián)合使用稱為Adams預(yù)測-校正系統(tǒng)，利用外插公式計算預(yù)測，用內(nèi)插公式進行校正。計算時需要注意以下兩點：1、 外插公式為顯式，內(nèi)插公式為隱式。故用內(nèi)插外插公式時需要進行迭代。2、 這種預(yù)測-校正法是四步法，計算Yn+1時，不但用到前一步信息，而且要用到更前三步信息，因此它不是自動開始的，實際計算時必須借助某種單步法，用Runge-Kutta格式為它提供初始值，依據(jù)局部截斷誤差公式得：進一步將Adams預(yù)測-校正系統(tǒng)加工成下列方案：運用上述方案計算時，要先一步的信息和

4、更前一步的信息因此在計算機之前必須給出初值和，可用其他單步法來計算，則一般令他為零。3. 源程序#include#include#include#include/微分方程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;/原函數(shù)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,dou

5、ble b,double h)int i;double c = (b-a)/h;double *y = new double c+1;for(i=0;i=c;i+)yi= 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;y0=y0;for(i=1;i=c;i+)yi= yi-1+ h* f(a+(i-1)*h,yi-1);/printf(%fn,f(a+(i-1)*h,yi-1);retu

6、rn 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;y0=y0;for(i=1;i=c;i+)yb=yi-1 +h* f(a+(i-1)*h,yi-1);yi = yi-1 + 0.5*h*( f(a+(i-1)*h,yi-1) + f(a+i*h,yb) );return y;/經(jīng)典K-R方法double* K_R(double a,double b,double h,double y0)double

7、 K1,K2,K3,K4,x;int i;double c = (b-a)/h ;double *y = new double c+1;y0=y0;for(i=1;i=c;i+)x=a+(i-1)*h;K1=f(x,yi-1);K2=f(x+0.5*h,yi-1+0.5*h*K1);K3=f(x+0.5*h,yi-1+0.5*h*K2);K4=f(x+h,yi-1+h*K3);yi=yi-1 + h*( K1 + 2*K2 + 2*K3 + K4)/6;return y;/4階Adams預(yù)測-校正方法double *Adams(double a,double b,double h,double

8、 y0)int i;double *y;double count = (b-a)/h;double dy100=0,x4=0;double p_0=0,p_1=0,c=0; double *r = new double *count+1;/動態(tài)初始化，儲存預(yù)測值和校值for(i=0;i=count;i+)ri = new double 2;if(r=NULL)printf(內(nèi)存分配失敗n);exit(0);for(i=0;icount;i+)ri0=0;ri1=0;y=K_R(a,b,h,y0);for(i=0;i4;i+)xi=a+h*i;dyi=f(xi,yi);ri0=yi;i=3;wh

9、ile(x3b)x3 = x3+h;p_1=y3+h*(55*dy3-59*dy2+37*dy1-9*dy0)/24; /預(yù)估c=p_1+251*(c-p_0)/270;/修正ri0=c;c=f(x3,c);/求fc=y3+h*(9*c+19*dy3-5*dy2+dy1)/24;/校正y3=c-19*(c-p_1)/270;/修正ri1 = y3;dy0=dy1;dy1=dy2;dy2=dy3;dy3=f(x3,y3);p_0=p_1;i+;return r;void main()int i , selet;double a=0, b=2, h=0.001, y0=1;while(1)prin

10、tf(請輸入下限， 上限， 步長 ， 初值：(用空格隔開)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、經(jīng)典K-R4、4階Adams預(yù)測-校正 5、修改數(shù)據(jù) 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(改進歐拉方法

11、：n);if(selet=3)y=K_R(a,b,h,y0);printf(經(jīng)典K-R方法(4階龍格-庫塔方法)：n);if(selet=4)double *r;r=Adams(a,b,h,y0);printf(4階Adams預(yù)測-校正方法:n);printf(x值 預(yù)估值 校正值 誤差n);printf(%0.3f %lf - -n,0,r00);printf(%0.3f %lf - -n,0.1,r10);printf(%0.3f %lf - -n,0.2,r20);for(i=3;i=10;i+)printf(%0.3f %f %f %0.20fn,a+i*h,ri0,ri1,fabs(ri0-ri1);if(selet=1|selet=2|selet=3)printf(x值 計算值 準(zhǔn)確值 誤差n);for(i=0;i改進歐拉歐拉方法4階adams預(yù)測校正方法2， 同一個方法的情況下，步長越小 結(jié)果準(zhǔn)確度越高；參考文獻：1 朱建新，李有法，數(shù)值計算方法(第三版)，高等教育出版社；2Adam預(yù)測-校正系統(tǒng)學(xué)生學(xué)習(xí)