微分方程-第一章習(xí)題5實(shí)驗(yàn)報(bào)告

1、第一章習(xí)題5實(shí)驗(yàn)報(bào)告一、實(shí)驗(yàn)?zāi)康耐ㄟ^(guò)用Euler法、中點(diǎn)格式、預(yù)報(bào)校正格式、Adams預(yù)報(bào)修正格式等算法求解一階常微分方程初值問(wèn)題的數(shù)值解，經(jīng)典四級(jí)四階R-K法是高精度單步法，四階Adams預(yù)估校正法是線性多步法，這兩個(gè)算法較實(shí)驗(yàn)一算法精度要高。通過(guò)本次實(shí)驗(yàn)，掌握R-K法的編程要領(lǐng)，掌握多步法初值的準(zhǔn)備方法，深刻領(lǐng)會(huì)微分方程數(shù)值解的實(shí)質(zhì)，體會(huì)單步法和線性多步法各自的優(yōu)缺點(diǎn)，熟練掌握各算法的計(jì)算機(jī)實(shí)現(xiàn)過(guò)程，并能從理論及實(shí)驗(yàn)結(jié)果分析各種算法的優(yōu)缺點(diǎn)。二、實(shí)驗(yàn)內(nèi)容用Euler法、中點(diǎn)格式、預(yù)報(bào)校正格式、經(jīng)典四級(jí)四階R-K格式、Adams預(yù)報(bào)修正格式算法求解一階常微分方程初值問(wèn)題的數(shù)值解，其精確解為。

2、三、算法1 Euler法其中，為步長(zhǎng)。2中點(diǎn)格式3預(yù)報(bào)校正格式4 四級(jí)四階格式：5 Adams預(yù)報(bào)修正格式其中為步長(zhǎng)。四 程序一、Euler法程序如下：#include <stdio.h>#define LEN 100int Euler(int a,int b,double h)double yLEN=0;double xLEN=0;int i=0;double num=0;num=(b-a)/h)+1;/節(jié)點(diǎn)個(gè)數(shù)/循環(huán)計(jì)算出各個(gè)節(jié)點(diǎn)的坐標(biāo)for(i=1;i<num;+i)x0=a;xi=a+i*h;/帶公式計(jì)算各個(gè)節(jié)點(diǎn)的y的值for(i=0;i<num;i+)y0=1

3、;yi+1=(yi/(2*xi)+(xi*xi)/(2*yi);/輸出x和y printf("n * Euler法求解方程 * n");printf("n xn ynnn");for(i=0;i<num;i+)printf(" %.4f %.4fn",xi,yi);printf("n * Euler法求解方程 * nn");return 0;二、中點(diǎn)格式程序如下：int Mid(int a,int b,double h)double yLEN=0;double xLEN=0;int i=0;double nu

4、m=0;num=(b-a)/h)+1;/節(jié)點(diǎn)個(gè)數(shù)/循環(huán)計(jì)算出各個(gè)節(jié)點(diǎn)的坐標(biāo)for(i=0;i<num;+i)x0=a;xi=a+i*h;/帶公式計(jì)算各個(gè)節(jié)點(diǎn)的y的值for(i=0;i<num;i+)y0=1;yi+1=yi+h*f(xi+h/2),(yi+(h/2)*(f(xi,yi);/輸出x和y printf("n * 中點(diǎn)格式求解方程 * n");printf(" xn ynnn");for(i=0;i<num;i+)printf(" %.4f %.4fn",xi,yi);printf("n * 中點(diǎn)

5、格式求解方程 * nn");return 0;三、預(yù)報(bào)校正格式程序如下：int PreCorr(int a,int b,double h)double y1LEN=0; double yLEN=0;double xLEN=0;int i=0;double num=0;num=(b-a)/h)+1;/節(jié)點(diǎn)個(gè)數(shù)/循環(huán)計(jì)算出各個(gè)節(jié)點(diǎn)的坐標(biāo)for(i=0;i<num;+i)x0=a;xi=a+i*h;/帶公式計(jì)算各個(gè)節(jié)點(diǎn)的y的值for(i=0;i<num;i+)y0=1;/初值為1 y1i+1=yi+h*f(xi,yi);yi+1=yi+(h/2)*(f(xi,yi)+f(xi+

6、1,y1i+1);/輸出x和y printf("n * 預(yù)報(bào)校正格式求解方程 * n");printf(" xn ynnn");for(i=0;i<num;i+)printf(" %.4f %.4fn",xi,yi);printf("n * 預(yù)報(bào)校正格式求解方程 * nn");return 0;四、經(jīng)典四級(jí)四階R-K程序如下：int RungeKutta(int a,int b,double h)double y1LEN=0; double yLEN=0;double xLEN=0;int i=0;doubl

7、e num=0;num=(b-a)/h)+1;/節(jié)點(diǎn)個(gè)數(shù)/循環(huán)計(jì)算出各個(gè)節(jié)點(diǎn)的坐標(biāo)for(i=0;i<num;+i)x0=a;xi=a+i*h;/帶公式計(jì)算各個(gè)節(jié)點(diǎn)的y的值for(i=0;i<num;i+)double K1=0;double K2=0;double K3=0;double K4=0;y0=1;/初值為1K1=f(xi,yi);K2=f(xi+h/2),(yi+h*K1/2);K3=f(xi+h/2),(yi+h*K2/2);K4=f(xi+h),(yi+h*K3);yi+1=yi+(h/6)*(K1+2*K2+2*K3+K4);/輸出x和y printf(&quo

8、t;n * 經(jīng)典四級(jí)四階R-K格式求解方程 * n");printf(" xn ynnn");for(i=0;i<num;i+)printf(" %.4f %.4fn",xi,yi);printf("n * 經(jīng)典四級(jí)四階R-K格式求解方程 * nn");return 0;六、Adams程序如下：int Adams(int a,int b,double h)double y1LEN = 0;double y2LEN = 0; double yLEN = 0;double xLEN = 0;int i = 0;double

9、 num = 0;double tmp = 0;num = (b-a)/h)+1;/節(jié)點(diǎn)個(gè)數(shù)/循環(huán)計(jì)算出各個(gè)節(jié)點(diǎn)的坐標(biāo)for(i = 0; i < num; +i)x0 = a;xi = a + i*h;for(i = 0; i < num; i+)y2i = sqrt(xi/2+(xi*xi)/2);y20 = 1;/帶公式計(jì)算各個(gè)節(jié)點(diǎn)的y的值for(i = 0; i < 4; i+)double K1 = 0;double K2 = 0;double K3 = 0;double K4 = 0;y0 = 1;/初值為1K1 = f(xi,yi);K2 = f(xi+h/2)

10、,(yi+h*K1/2);K3 = f(xi+h/2),(yi+h*K2/2);K4 = f(xi+h),(yi+h*K3);yi+1 = yi+(h/6)*(K1+2*K2+2*K3+K4);y1i = yi;for(i = 4;i < num;i+) y1i+1 = yi+(h/24)*(55*f(xi,yi)-59*f(xi-1,yi-1)+37*f(xi-2,yi-2)-9*f(xi-3,yi-3); yi+1 = yi+(h/24)*(9*f(xi+1,y1i+1)+19*f(xi,yi)-5*f(xi-1,yi-1)+f(xi-2,yi-2);/輸出x和y printf(&q

11、uot;n *Adams預(yù)報(bào)校正格式求解方程 * n");printf(" xn yn y(精確解)nn");for(i = 0; i < num; i+)printf(" %.4f %.4f %.4f n",xi,yi,y2i);printf("n *Adams預(yù)報(bào)校正格式求解方程 * nn");return 0;五、精確解程序如下：int ac_solve(int a,int b,double h)double yLEN=0;double xLEN=0;int i=0;double num=0;num=(b-a)/

12、h)+1;/節(jié)點(diǎn)個(gè)數(shù)/循環(huán)計(jì)算出各個(gè)節(jié)點(diǎn)的坐標(biāo)for(i=1;i<num;+i)x0=a;xi=a+i*h;/帶公式計(jì)算各個(gè)節(jié)點(diǎn)的y的值for(i=0;i<num;i+)yi=sqrt(xi/2+(xi*xi)/2);y0=1;/輸出x和y printf("n * 精確解 *n");printf("n xn ynnn");for(i=0;i<num;i+)printf(" %.4f %.4fn",xi,yi);printf("n * 精確解 *nn");return 0;調(diào)用前四個(gè)函數(shù)#inclu

13、de<stdio.h>#include<math.h>#include "Euler.h"#define LEN 100int main()int a=1;/區(qū)間左端點(diǎn)int b=2;/區(qū)間右端點(diǎn)int case_flag=0;double h=0.1;/步長(zhǎng)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);/

14、中點(diǎn)格式break;case 3:PreCorr(a,b,h);/預(yù)報(bào)校正格式break;case 4:RungeKutta(a,b,h);/經(jīng)典四級(jí)四階R-K格式break;case 5:ac_solve(a,b,h);/精確解break;default:printf("沒(méi)有您所選擇的序號(hào)!n");break;return 0;調(diào)用Adams預(yù)報(bào)校正格式主函數(shù)int main(void)int a = 1;/區(qū)間左端點(diǎn)int b = 2;/區(qū)間右端點(diǎn)double h = 0.1;/步長(zhǎng)Adams(a,b,h);/Adams法return 0;五、數(shù)值結(jié)果及分析運(yùn)算結(jié)果可用

15、下面表格表示： 表5.1四種方法以及方程精確解的計(jì)算結(jié)果表xiEuler格式中點(diǎn)格式預(yù)報(bào)-校正格式經(jīng)典四級(jí)四階R-K格式Adams格式精確解1.0000 1.0000 1.0000 1.0000 1.0000 1.00001.0000 1.1000 1.0000 1.10251.10251.10251.10251.07471.2000 1.05951.2100 1.2100 1.2100 1.21001.14891.3000 1.1210 1.32231.32231.3223 1.32231.22271.4000 1.18491.43951.43951.4394 1.43941.29611.5

16、000 1.25021.56131.56131.5612 1.56121.36931.6000 1.31661.68771.68761.6876 1.68761.44221.7000 1.38361.81841.81841.8184 1.81841.51491.8000 1.45131.95351.95351.9535 1.95351.58751.9000 1.51942.09282.09282.0927 2.09271.65982.0000 1.58782.23622.23612.2361 2.23611.7321 表5.2四種方法和方程精確解的差值表xiEuler格式中點(diǎn)格式預(yù)報(bào)-校正格式

17、經(jīng)典四級(jí)四階R-K格式Adams格式精確解1.0000 0.0000 0.0000 0.0000 0.0000 0.0000 1.0000 1.1000 0.0747-0.0278-0.0278-0.0278-0.02781.07471.2000 0.0894-0.0611 -0.0611 -0.0611 -0.0611 1.14891.3000 0.1017 -0.0996-0.0996-0.0996-0.09961.22271.4000 0.1112-0.1434-0.1434-0.1433-0.14331.29611.5000 0.1191-0.192-0.192-0.1919-0.19

18、191.36931.6000 0.1256-0.2455-0.2454-0.2454-0.24541.44221.7000 0.1313-0.3035-0.3035-0.3035-0.30351.51491.8000 0.1362-0.366-0.366-0.366-0.3661.58751.9000 0.1404-0.433-0.433-0.4329-0.43291.65982.0000 0.1443-0.5041-0.504-0.504-0.5041.7321求和1.1739 -2.3760 -2.3758 -2.3755 -2.3755 由以上表格可知：從表5.1和表5.2可以看出在四種

19、方法中Euler格式計(jì)算的方程的解與精確解相差最大，即誤差較大；中點(diǎn)格式，預(yù)報(bào)-校正格式和經(jīng)典四級(jí)四階R-K格式的誤差較小，其中經(jīng)典四級(jí)四階R-K格式離精確解最近，所以，經(jīng)典四級(jí)四階R-K格式方法在前四種方法中求解方程的解的誤差最小，效果最好。通過(guò)計(jì)算結(jié)果，可知預(yù)報(bào)格式已是很好的近似，同時(shí)Adams內(nèi)插法又有較高的精度，故我們用相應(yīng)的外插公式作為預(yù)報(bào)格式，然后用內(nèi)插格式進(jìn)行迭代，從而形成Adams預(yù)報(bào)修正格式，所以Adams預(yù)報(bào)修正格式效果顯著，能達(dá)到較高的精度，盡管只迭代一次卻與真解十分相近。第一章習(xí)題9實(shí)驗(yàn)報(bào)告一、實(shí)驗(yàn)?zāi)康耐ㄟ^(guò)用Euler法求解高階常微分方程初值問(wèn)題的數(shù)值解，深刻領(lǐng)會(huì)微分方

20、程數(shù)值解的實(shí)質(zhì)，熟練掌握該算法的計(jì)算機(jī)實(shí)現(xiàn)過(guò)程，并能從理論及實(shí)驗(yàn)結(jié)果分析該算法的優(yōu)缺點(diǎn)。 通過(guò)本次實(shí)驗(yàn)，學(xué)會(huì)將一階常微分方程初值問(wèn)題的數(shù)值解法推廣到一階常微分方程組初值問(wèn)題中。這里以Euler法為例進(jìn)行推廣，體會(huì)其相同之處及不同之處，同時(shí)學(xué)會(huì)向量值函數(shù)的表達(dá)。二、實(shí)驗(yàn)內(nèi)容用Euler法求解高階常微分方程初值問(wèn)題的數(shù)值解，其精確解為。三、算法 Euler法求解高階常微分方程初值問(wèn)題其中為步長(zhǎng)。四、程序4.1頭文件 功能： 頭文件中包含了四個(gè)方法以及方程精確解求解的具體實(shí)現(xiàn)過(guò)程，并將各個(gè)方法打包成函數(shù)，方便在主文件中進(jìn)行調(diào)用。各個(gè)函數(shù)分別為double f1(double x,double y1,

21、double y2);double f2(double x,double y1,double y2);頭文件程序：double f1(double x,double y1,double y2)return y2;double f2(double x,double y1,double y2)return -y1;4.2主文件主文件程序：int main(void)double xARRARY_LEN = 0;double yARRARY_LENARRARY_LEN = 0;double y1ARRARY_LENARRARY_LEN = 0;int i=0;int a = 0;int b = 1;

22、double h = 0.1;double num=0;num=(b-a)/h)+1;/節(jié)點(diǎn)個(gè)數(shù)/循環(huán)計(jì)算出各個(gè)節(jié)點(diǎn)的坐標(biāo)for(i=1;i<num;+i)x0=a;xi=a+i*h;/計(jì)算精確解for(i=0;i<num;i+)y100=1;y110=1;y10i+1=sin(xi+1);y11i+1=cos(xi+1);/帶公式計(jì)算各個(gè)節(jié)點(diǎn)的y的值for(i=0;i<num;i+)y00=0;y10=1;y0i+1=y0i+h*f1(xi,y0i,y1i);y1i+1=y1i+h*f2(xi,y0i,y1i);/輸出x和y printf("n * Euler法求解高階方程 * n");printf("n xn y1n