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

數(shù)值分析常微分?jǐn)?shù)值解本科實(shí)習(xí)課程學(xué)院名稱專業(yè)名稱學(xué)生姓名學(xué)生學(xué)號(hào)指導(dǎo)教師實(shí)驗(yàn)地點(diǎn)實(shí)驗(yàn)成績生課程設(shè)計(jì)報(bào)告數(shù)值分析管理科學(xué)學(xué)院信息與計(jì)算科學(xué)二〇一六年六月二〇一六年六常微分方程數(shù)值解法是計(jì)算數(shù)學(xué)的一個(gè)分支.是解常微分方程各類定解問題的數(shù)值方法.現(xiàn)有的解析方法只能用于求解一些特殊類型的定解問題，實(shí)用上許多很有價(jià)值的常微分方程的解不能用初等函數(shù)來表示，常常需要求其數(shù)值解.所謂數(shù)值解，是指在求解區(qū)間內(nèi)一系列離散點(diǎn)處給出真解的近似值.這就促成了數(shù)值方法的產(chǎn)生與發(fā)展.關(guān)鍵詞：數(shù)值解法；常微分1.實(shí)驗(yàn)內(nèi)容和要求常微分方程初值問題式列出10個(gè)等距節(jié)點(diǎn)上的計(jì)算值和精確值，并比較他們的計(jì)算精確度。其中多步法需要的初值由經(jīng)典R-K法提供。運(yùn)算時(shí)，取足以表示計(jì)算精2.算法說明函數(shù)及變量說明函數(shù)/變量名fp_fEuler，Euler_ProK_R,AdamsMain類型doubledoubledouble函數(shù)及變量定義說明要求的微分方程微分方程原函數(shù)msVoid主Void4i4ii3y(x)=y(x)+hf(x,y(x))i+1iii(i=0，1，2，3，......n-(x)=y(x)+hf(x,y(x))i+1iiiy(x)=y(x)+h[f(x,y(x))+f(x,(x))]i+1i2iii+1i+1(i=0，1，2......n-1)(其yi+1=y()|3=i+2i+22|K=f(x+h,y+hK)Adsms式聯(lián)合使用稱為Adams預(yù)測-校正系統(tǒng)，利用外插公式計(jì)算預(yù)測，用內(nèi)插公式進(jìn)行校正。計(jì)算時(shí)需要注意以下兩點(diǎn)：1、外插公式為顯式，內(nèi)插公式為隱式。故用內(nèi)插外插公式時(shí)需要進(jìn)2、這種預(yù)測-校正法是四步法，計(jì)算Yn+1時(shí)，不但用到前一步信息，而且要用到更前三步信息f，ff，因此它不是自動(dòng)開始的，實(shí)nnn-3際計(jì)算時(shí)必須借助某種單步法，用Runge-Kutta格式為它提供初始值y,y,y，依據(jù)局部截?cái)嗾`差公式得：運(yùn)用上述方案計(jì)算y時(shí)，要先一步的信息y,y',pc和更前一步的信息n+1nnnnynn因此在計(jì)算機(jī)之前必須給出初值y和p-c，y可用其他單步法來計(jì)算，p-c則一般令他為零。3.源程序#include<>#include<>#include<>#include<>//微分方程doublef(doublex,doubley){return(-y+cos(2*x)-2*sin(2*x)+2*x*exp(-x));//returnx*pow(y,-2)*;}//原函數(shù)doublep_f(doublex){returnx*x*exp(-x)+cos(2*x);//returnpow+pow(x,2),;}//求精確解double*Exact(doublea,doubleb,doubleh){inti;doublec=(b-a)/h;double*y=newdouble[c+1];for(i=0;i<=c;i++)y[i]=p_f(a+i*h);returny;}//歐拉方法double*Euler(doublea,doubleb,doubleh,doubley0){inti;doublec=(b-a)/h;double*y=newdouble[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]));}returny;}//改進(jìn)歐拉方法double*Euler_Pro(doublea,doubleb,doubleh,doubley0){inti;doublec=(b-a)/h,yb;double*y=newdouble[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]+*h*(f(a+(i-1)*h,y[i-1])+f(a+i*h,yb));}returny;}double*K_R(doublea,doubleb,doubleh,doubley0){doubleK1,K2,K3,K4,x;inti;doublec=(b-a)/h;double*y=newdouble[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+*h,y[i-1]+*h*K1);K3=f(x+*h,y[i-1]+*h*K2);K4=f(x+h,y[i-1]+h*K3);y[i]=y[i-1]+h*(K1+2*K2+2*K3+K4)/6;}returny;}double**Adams(doublea,doubleb,doubleh,doubley0){inti;double*y;doublecount=(b-a)/h;doubledy[100]={0},x[4]={0};doublep_0=0,p_1=0,c=0;double**r=newdouble*[count+1];//動(dòng)態(tài)初始化，儲(chǔ)存預(yù)測值和校值for(i=0;i<=count;i++)r[i]=newdouble[2];if(r==NULL){printf("內(nèi)存分配失敗\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];}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;//預(yù)估c=p_1+251*(c-p_0)/270;//修正r[i][0]=c;c=f(x[3],c);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;}returnr;}voidmain(){inti,selet;doublea=0,b=2,h=,y0=1;while(1){開)\n");scanf("%lf%lf%lf%lf",&a,&b,&h,&y0);doublec=(b-a)/h;double*yx,*y;yx=Exact(a,b,h);while(1){Adams預(yù)測-校正5、修改數(shù)據(jù)6、退出\n");scanf("%d",&selet);system("cls");if(selet==1){y=Euler(a,b,h,y0);}if(selet==2){y=Euler_Pro(a,b,h,y0);}if(selet==3){y=K_R(a,b,h,y0);\n");}if(selet==4){double**r;r=Adams(a,b,h,y0);printf("x值預(yù)估值校正值誤差\n");printf("%%lf--\n",0,r[0][0]);printf("%%lf--\n",,r[1][0]);printf("%%lf--\n",,r[2][0]);for(i=3;i<=10;i++)printf("%%f%f%\n",a+i*h,r[i][0],r[i][1],fabs(r[i][0]-r[i][1]));}{printf("x值計(jì)算值準(zhǔn)確值誤差\n");for(i=0;i<=10;i++)printf("%-y[i]));}}}%lf%lf%\n",(a+i)*h,y[i],yx[i],fabs(yx[i]}if(selet==5)break;if(selet==6)exit(1);continue;4.實(shí)驗(yàn)結(jié)果1、改進(jìn)歐拉方法長精確解長精確解長精精確解5.對(duì)比分析改進(jìn)歐拉方法改進(jìn)歐拉方法經(jīng)典K-R方法4階adams預(yù)測-校正方法歐拉方法方法誤步差長1，在步長相同的時(shí)候，結(jié)果準(zhǔn)確度經(jīng)典K-R>