常微分实习报告

1、课程实习报告课程名称：常微分方程课程实习实习题目：常微分方程数值求解问题的实习姓 名:系: 专 业: 年 级: 学 号: 指导教师： 职 称:课程实习报告结果评定评语：成绩：指导教师签字：评定日期：1.实习的目的和任务 1.2.实习要求 1.3.实习地点 1.4.主要仪器设备 1.5.实习内容 1.5.1用不同格式对同一个初值问题的数值求解及其分析 .15.1.1 求精确解1.5.1.2用欧拉法求解6.5.1.3用改进欧拉法求解 9.5.1.4用4级4阶龙格库塔法求解 1.25.1.5 问题讨论与分析 155.2 一个算法不同不长求解同一个初值问题及其分析 .186.结束语 2.9.参考文献

2、2.9.常微分方程课程实习1.实习的目的和任务目的：通过课程实习能够应用MATLAB来计算微分方程(组)的数值解；了解常微分方程数值解。任务：通过具体的问题，利用MATLAB件来计算问题的结果，分析问题的结论。2. 实习要求能够从案例的自然语言描述中，抽象出其中的数学模型；能够熟练应用 所学的数值解计算方法；能够熟练使用MATLAB软件；对常微分方程数值解有所认识，包括对不同算法有所认识和对步长有所认识。3. 实习地点数学实验室、学生宿舍4. 主要仪器设备计算机、Microsoft Win dows 7Matlab 7.05. 实习内容5.1用欧拉方法，改进欧拉方法，4阶龙格一库塔方法分别求下

3、面微分方程的初值 dy/dx=y*cos(x+2) y(-2)=1 x -2，05.1.1求精确解 变量分离方程情形：形如3 = f (x)* g(y)的方程，这里f (x), g(y)分别是x, y的dx连续函数.如果g(y),我们可将方程改写成 型 f(x)dx,这样,变量就”分g(y)离”开来了 ,两边同时积分即可：-业f(x)dx c,c为任意常数.g(y) 常数变易法：一阶线性微分方程 鱼=f (x)* y g(x),其中f (x), g(x)在考虑区dx间上是的连续函数.可先解出方程 鱼二f(x)* y的解,这是属于变量分dx离方程情形，可解得：y =c*exp( f (x)dx)

4、,这里c是任意常数.然后将变c得：dy dc(x)*exp( f(x)dx) c(x)* f(x)*exp( f(x)dx) = f (x)* c(x)*exp( f (x)dx) g(x)dx dx所以可解得 c(x)二 f(x)*exp(- f (x)dx)dx* cl,这里 cl 是任意常数.将 c(x) = j f (x)*exp(- f (x)dx)dx - cl 代入y =c(x)*exp( f (x)dx)可得原方程的通解：y =( f(x)*exp(- f(x)dx)dx c1)*exp( f (x)dx),cl 为任意常数.恰当微分方程情形:形如m(x, y)dx - n(x

5、, y)dy二0的一阶微分方程,这里 假设m(x, y), n(x, y)在某矩形域内是的连续函数,且具有连续的一阶偏导数 若 卫=口，则为恰当微分方程.判断为恰当微分方程后,则可用如下解法:设u(x, y)=c是原方程的解，则 =m，所以设m(x, y)dx v(y)，；：( .m(x,y)dx v(y)_ ；:( m(x, y)dx) dv(y)=，所以dy一(rn(x, y)dx) =dv(y)，由此 v(y)二n m(x,y)dxdy ,由此可解得dym(x, y)dx 亠 i n -一 m(x, y)dxdy = c, c 为任意常数首先可以求得其精确解为：y=exp(s in (x

6、+2) x=-2:0.1:2; y= exp(si n( x+2) plot(x,y,r.-);Columns 1 through 41.000000000000001.34382524373165Columns 5 through 81.476121946445731.904496534386731.104986830331691.615146296442081.219778556000621.75881884576699Columns 9 through 122.049008650164272.438071505156332.188741912604612.31977682471585Co

7、lumns 13 through 162.539682532380782.711481017682162.621005926286702.67901644757271Columns 17 through 202.717123008431282.695718599203822.648113847390782.57616043684702Columns 21 through 242.482577728015002.10792744704554Columns 25 through 281.964942888556771.53323499677732Columns 29 through 321.397

8、923819945001.04245724559826Columns 33 through 360.943296953050410.70413637458179Columns 37 through 400.642415207750370.502697762537372.370757126170311.819336991081061.270295218811560.854066948581020.58870142585634Column 410.46916418587400Data =2.244530580577551.674477827371161.151562836514530.774497

9、302497390.54234231959371-1.80000000000000-1.70000000000000-1.60000000000000-1.500000000000001.219778556000621.343825243731651.476121946445731.61514629644208-1.40000000000000-1.300000000000001.758818845766991.904496534386731.200000000000002.049008650164271.100000000000002.18874191260461-1.00000000000

10、000 2.31977682471585-0.90000000000000 2.43807150515633-0.800000000000002.53968253238078-0.70000000000000 2.62100592628670-0.60000000000000 2.67901644757271 -0.50000000000000 2.71148101768216-0.40000000000000 2.71712300843128 -0.30000000000000 2.69571859920382-0.20000000000000 2.64811384739078-0.1000

11、00000000002.576160436847020 2.482577728015000.10000000000000 2.370757126170310.20000000000000 2.244530580577550.30000000000000 2.107927447045540.400000000000000.500000000000000.600000000000000.700000000000000.800000000000000.900000000000001.964942888556771.819336991081061.674477827371161.53323499677

12、7321.397923819945001.100000000000001.200000000000001.300000000000001.400000000000001.500000000000001.600000000000001.700000000000001.800000000000001.900000000000002.000000000000001.042457245598260.943296953050410.854066948581020.774497302497390.704136374581790.642415207750370.588701425856340.5423423

13、19593710.502697762537370.469164185874000.52.521.51亠0-2-1.5-1-0.50.51.55.1.2用欧拉法求解设常微分方程的初始问题史=f(x, y),a - x -b,(1)dxy(a) = .(2)有唯一解。则由欧拉法求初值问题（1）, （2）的数值解的差分方程为:yn 1 二 yn hf（Xn,yn）,yi 二.程序如下：建立函数文fl.mfun cti on x,y=f1(fu n, x_spa n,yO,h)x=x_spa n(1):h:x_spa n(2);y(i)=y0;for n=1:le ngth(x)-1y( n+1)=y

14、( n)+h*feval(fu n,x( n),y( n);endx=x;y=y;在MATLAB输入以下程序： clear all fun=inlin e( y*cos(x+2); x,y1=f1(fu n,-2,2,1,0.1); x,y1 plot(x,y1,g*-)结果及其图象：ans =-2.00000000000000-1.90000000000000-1.80000000000000-1.70000000000000-1.60000000000000-1.50000000000000-1.400000000000001.000000000000001.100000000000001

15、.209450458180581.327984655342341.454851875167081.588852606593921.728287540690011.100000000000001.000000000000000.900000000000000.800000000000000.700000000000000.600000000000000.500000000000000.400000000000000.300000000000000.200000000000000.1000000000000000.100000000000000.200000000000000.3000000000

16、00000.400000000000000.500000000000000.600000000000000.700000000000000.800000000000000.900000000000001.000000000000001.100000000000001.200000000000001.300000000000001.400000000000001.500000000000001.600000000000002.154344360817052.288260553794212.411895799158422.521298457136512.612659661865862.682548

17、001780242.728142503735772.747440620382272.739418225015642.704122329429032.642684103673772.557248883750392.450829780426752.327100593658172.190150463724832.044225990027361.893486050208131.741790624182991.592538544524401.448561571205111.312074863782761.184677883555021.067395661994220.960748409479460.86

18、4837397675810.779436454229260.704080678711320.638146572714180.580920241720551.900000000000000.489600406402422.000000000000000.4540587312867232.521.510.50-2-1.5-1-0.500.511.525.1.3用改进欧拉法求解：计算公式为：y0）1 *n hf （Xn,yn）,=% pf（Xn,yn） + f （人十朋）,y = ,n= 0,1,2,.,k = 0,1,2,.,即先用欧拉法得（xn,yn）,进而由（3）的第一式得初始近似值yn0）1

19、 ,然后再用（3） 的第二式进行迭代，反复改进这个近似值，直到| ynT - yk ：；（：为所允许的误差）为止，并把ynk1取作为y（xn 1）的近似值ynv，这个方法就称为改进欧拉法。 通常称（3）为预报校正公式，其中第一式称为预报公式，第二式称为校正公式。这个公式还可以写为：（下文改进欧拉法计算就以下面为准）yn 厂 yn2(klk2),ki =f (Xnn),k2 =f 区 h, ynhki),y。=.程序如下：建立函数文件f2.mfun cti on x,y=f2(fu n, x_spa n,yO,h)x=x_spa n(1):h:x_spa n(2);y(1)=yO;for n=1

20、:le ngth(x)-1k1=feval(fu n,x( n),y( n);y( n+1)=y( n)+h*k1; k2=feval(fu n,x(n+1),y( n+1);y( n+1)=y( n)+h*(k1+k2)/2; end x=x;y=y;在MATLAB输入以下程序： clear all fun=i nlin e( y*cos(x+2); x,y2=f2(fu n,-2,2,1,O.1); x,y2 plot(x,y2,b+-)结果及其图象：ans =-2.00000000000000-1.90000000000000-1.800000000000001.000000000000

21、001.104725229090291.219207229363991.500000000000001.400000000000001.300000000000001.200000000000001.100000000000001.000000000000000.900000000000000.800000000000000.700000000000000.600000000000000.500000000000000.400000000000000.300000000000000.200000000000000.1000000000000000.100000000000000.2000000

22、00000000.300000000000000.400000000000000.500000000000000.600000000000000.700000000000000.800000000000000.900000000000001.000000000000001.100000000000001.200000000000001.613388617478771.756604945660731.901814952758942.045861837217642.185146633762432.315763555724452.433682968959042.534971671734682.616

23、033679010132.673849667845002.706190767907232.711783264044682.690405219408782.642903530441302.571129346284082.477799560853782.366300571501322.240456332742542.104285124930581.961768315515181.816650280351961.672282601880951.531518885262901.396660163895361.269445745893071.151080940502511.042291474098520

24、.943394336221440.854375884572411.500000000000001.600000000000001.700000000000001.800000000000001.900000000000002.000000000000000.704729719126620.643092734337550.589432936264430.543103926871230.503471430640500.469937065943502.5+十4-十七 专 +1.54-七+卡0.50 L-2-1.5-1-0.50.51.55.1.4用4阶龙格一库塔求解标准的四阶龙格-库塔公式(亦称为经

25、典的四阶龙格-库塔公式)yn 1 =yn 一(匕 2k22& kQ6 2kf (Xn,yn),丄1丄hk2 =f (Xn ：h, yn -k1),221 hk3 = f (Xn ：h, yn -k2),2 2k f (Xnh,ynhka),公式的截断误差阶为 o(h5 )程序如下：建立函数文件 f3.mfunction x,y=f3(fun,x_span,y0,h)x=x_span (1):h:x_span(2);y(1)=y0;for n=1:length(x)-1k1=feval(fun,x(n),y(n); k2=feval(fun,x(n)+h/2,y(n)+h/2*k1); k3=f

26、eval(fun,x(n)+h/2,y(n)+h/2*k2);k4=feval(fun,x(n+1),y(n)+h*k3); y(n+1)=y(n)+h*(k1+2*k2+2*k3+k4)/6;endx=x;y=y;在 MATLAB 输入以下程序： clear all; fun=inline( y*cos(x+2); x,y3=f3(fun,-2,2,1,0.1); x,y3 plot(x,y3, y*-) 结果及其图象： ans =-2.00000000000000-1.90000000000000-1.80000000000000-1.70000000000000-1.6000000000

27、00001.000000000000001.104986745696811.219778373518281.343824953625391.400000000000001.300000000000001.200000000000001.100000000000001.000000000000000.900000000000000.800000000000000.700000000000000.600000000000000.500000000000000.400000000000000.300000000000000.200000000000000.1000000000000000.10000

28、0000000000.200000000000000.300000000000000.400000000000000.500000000000000.600000000000000.700000000000000.800000000000000.900000000000001.000000000000001.100000000000001.200000000000001.300000000000001.758818219617671.904495806501982.049007830855662.188741013447442.319775857524332.438070481408852.5

29、39681463116192.621004822310582.679015319711092.711479876846582.717121865404832.695717464250942.648112729936422.576159345482252.482576670951542.370756112042032.244529619279342.107926550209661.964942069343181.819336263174131.674477203345411.533234486212941.397923427764521.270294944247731.1515626729423

30、11.042457181239250.943296972359340.854067034002240.774497436252561.600000000000001.700000000000001.800000000000001.900000000000002.000000000000000.642415391183080.588701616051610.542342508643010.502697945435360.469164360045032.51.50.50 L-2rrr-1.5-1-0.500.511.55.1.5问题讨论与分析由以上数值分析结果绘制表格:精确解欧拉万法改进的欧拉方法

31、四阶龙格-库塔方法xiyiyi误差yi误差yi误差-2.001.0000001.0000000.0000001.0000000.0000001.0000000.000000-1.901.1049871.1000000.0049871.1047250.0002621.1049870.000000-1.801.2197791.2094500.0103281.2192070.0005711.2197780.000000-1.701.3438251.3279850.0158411.3428980.0009271.3438250.000000-1.601.4761221.4548520.0212701.

32、4747970.0013251.4761220.000000-1.501.6151461.5888530.0262941.6133890.0017581.6151460.000001-1.401.7588191.7282880.0305311.7566050.0022141.7588180.000001-1.301.9044971.8709290.0335671.9018150.0026821.9044960.000001-1.202.0490092.0140260.0349832.0458620.0031472.0490080.000001-1.102.1887422.1543440.034

33、3982.1851470.0035952.1887410.0000011.002.3197772.2882610.0315162.3157640.0040132.3197760.0000010.902.4380722.4118960.0261762.4336830.0043892.4380700.0000010.802.5396832.5212980.0183842.5349720.0047112.5396810.0000010.702.6210062.6126600.0083462.6160340.0049722.6210050.0000010.602.6790162.6825480.003

34、5322.6738500.0051672.6790150.0000010.502.7114812.7281430.0166612.7061910.0052902.7114800.0000010.402.7171232.7474410.0303182.7117830.0053402.7171220.0000010.302.6957192.7394180.0437002.6904050.0053132.6957170.0000010.202.6481142.7041220.0560082.6429040.0052102.6481130.0000010.102.5761602.6426840.066

35、5242.5711290.0050312.5761590.0000010.002.4825782.5572490.0746712.4778000.0047782.4825770.0000010.102.3707572.4508300.0800732.3663010.0044572.3707560.0000010.202.2445312.3271010.0825702.2404560.0040742.2445300.0000010.302.1079272.1901500.0822232.1042850.0036422.1079270.0000010.401.9649432.0442260.079

36、2831.9617680.0031751.9649420.0000010.501.8193371.8934860.0741491.8166500.0026871.8193360.0000010.601.6744781.7417910.0673131.6722830.0021951.6744770.0000010.701.5332351.5925390.0593041.5315190.0017161.5332340.0000010.801.3979241.4485620.0506381.3966600.0012641.3979230.0000000.901.2702951.3120750.041

37、7801.2694460.0008491.2702950.0000001.001.1515631.1846780.0331151.1510810.0004821.1515630.0000001.101.0424571.0673960.0249381.0422910.0001661.0424570.0000001.200.9432970.9607480.0174510.9433940.0000970.9432970.0000001.300.8540670.8648370.0107700.8543760.0003090.8540670.0000001.400.7744970.7794360.004

38、9390.7749700.0004730.7744970.0000001.500.7041360.7040810.0000560.7047300.0005930.7041370.0000001.600.6424150.6381470.0042690.6430930.0006780.6424150.0000001.700.5887010.5809200.0077810.5894330.0007320.5887020.0000001.800.5423420.5316520.0106900.5431040.0007620.5423430.0000001.900.5026970.4896000.013

39、0970.5034710.0007740.5026970.0000002.000.4691640.4540580.0151060.4699370.0007730.4691640.000000 plot(x,y,r+-) hold on,plot(x,y1,b-) plot(x,y,r+-) hold on,plot(x,y2,b-) plot(x,y,r+-) hold on,plot(x,y3,b-)精确解与欧拉方程比较2.51.50.5十十+十%4-4-0 1L00.51.5-2-1.5-1-0.5精确解与改进后的欧拉方程比较3由上表和图可以看出欧拉法误差最大，而改进欧拉和龙格一库塔方法误

40、差相 对较小，并且龙格一库塔方法误差最小且大部分值都跟精确值相同。由欧拉图与精确图相比可清晰看到，随着X的增加，函数值与精确值的偏差越来越大。5.2选择用欧拉方法，改进欧拉方法，4阶龙格一库塔方法之一取不同步长分别 求下面微分方程的初值dy/dx=y*(cosx+1) y(-2)=1 x -2，0，并对结果进行 分析说明，给出你的结论。使用欧拉方法：取步长h=0.05,如下：程序如下：建立函数文f4.mfun cti on x,y=f4(fu n, x_spa n,yO,h)x=x_spa n(1):h:x_spa n(2);y(1)=y0;for n=1:le ngth(x)-1y(n+1)

41、=y( n)+h*feval(fu n, x( n),y( n);endx=x;y=y;在 MATLAB 输入以下程序： clear all fun=inline( y*cos(x+2) ); x,y4=f4(fun,-2,0,1,0.05); x,y4 plot(x,y4,g+-)ans =-2.00001.0000-1.95001.0500-1.90001.1024-1.85001.1573-1.80001.2145-1.75001.2740-1.70001.3357-1.65001.3995-1.60001.4653-1.55001.5327-1.50001.6018-1.45001.6

42、720-1.40001.7433-1.35001.8153-1.30001.8875-1.25001.9597-1.20002.0314-1.15002.1021-1.10002.1715-1.05002.2390-0.9500 2.3664-0.9000 2.4252-0.8500 2.4803-0.8000 2.5309-0.7500 2.5768-0.7000 2.6174-0.6500 2.6524-0.6000 2.6814-0.5500 2.7042-0.5000 2.7205-0.4500 2.7301-0.4000 2.7330-0.3500 2.7290-0.3000 2.7

43、182-0.2500 2.7007-0.2000 2.6766-0.1500 2.6462-0.1000 2.6097-0.0500 2.56760 2.52002.62 42.2re1.61412121.8-1.6.4-1.2-1-0.8-oe-0.4-0.2取步长h=0.02程序如下：建立函数文f5.mfun cti on x,y=f5(fu n, x_spa n,y0,h)x=x_spa n(1):h:x_spa n(2);y(i)=yo；for n=1:le ngth(x)-1y( n+1)=y( n)+h*feval(fu n,x( n),y( n);endx=x;y=y;在MATL

44、AB输入以下程序： clear all fun=inlin e( y*cos(x+2); x,y5=f5(fu n,-2,0,1,0.02); plot(x,y5,叶-)-1.9800-1.9600-1.9400-1.9200-1.9000-1.8800-1.8600-1.8400-1.8200-1.8000-1.7800-1.7600-1.7400-1.7200-1.7000-1.6800-1.6600-1.6400-1.6200-1.6000-1.5800-1.5600-1.5400-1.5200-1.5000-1.4800-1.46001.02001.04041.06121.08241.

45、10401.12591.14831.17101.19411.21761.24151.26571.29031.31531.34051.36621.39211.41831.44491.47171.49881.52621.55381.58171.60971.63801.6664-1.4000-1.3800-1.3600-1.3400-1.3200-1.3000-1.2800-1.2600-1.2400-1.2200-1.2000-1.1800-1.1600-1.1400-1.1200-1.1000-1.0800-1.0600-1.0400-1.0200-1.0000-0.9800-0.9600-0.

46、9400-0.9200-0.9000-0.8800-0.86001.75261.78151.81051.83951.86861.89771.92671.95561.98452.01332.04192.07042.09862.12662.15442.18182.20902.23572.26212.28812.31352.33852.36302.38692.41032.43302.45512.4765-0.8000-0.7800-0.7600-0.7400-0.7200-0.7000-0.6800-0.6600-0.6400-0.6200-0.6000-0.5800-0.5600-0.5400-0

47、.5200-0.5000-0.4800-0.4600-0.4400-0.4200-0.4000-0.3800-0.3600-0.3400-0.3200-0.3000-0.2800-0.26002.53632.55472.57222.58892.60482.61972.63372.64682.65892.67002.68012.68932.69732.70442.71042.71532.71912.72192.72352.72412.72362.72202.71942.71562.71082.70492.69792.6899-0.22002.6707-0.20002.6596-0.18002.6

48、475-0.16002.6345-0.14002.6205-0.12002.6055-0.10002.5897-0.08002.5729-0.06002.5553-0.04002.5369-0.02002.517602.4976以下为比较h=0.1、0.05、0.02在取相同x值，比较初值:精确值H=0.1H=0.05H=0.02xiyiyi误差yi误差yi误差-2.01.00001.000001.000001.00000-1.91.10501.10000.00501.10240.00261.10400.0010-1.81.21981.20950.01031.21450.00531.21760.0022-1.71.34381.32800.01581.33570.00711.34050.0035-1.61.47611.45490.0