常微分实习报告-常微分方程数值求解问题

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

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

4、( J f (x)dx),这里c是任意常数.然后将变c易为x的待定函数c(x),令y =c(x)*exp( j f (x)dx),将其代入原方程可得：二 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* c1,这里 cl 是任意常数.将 c(x) = J f (x)*exp( - J f (x)dx)dx + c1 代入y =c(x)*exp( jf (x)dx)可得原方程的通解：y =( J f (x

5、)*exp( f (x)dx)dx+c1)*exp( J f (x)dx), cl 为任意常数.恰当微分方程情形：形如m(x, y)dx+n(x, y)dy = 0的一阶微分方程，这里假设m(x, y), n(x, y)在某矩形域内是的连续函数，且具有连续的一阶偏导数.若如=2，则为恰当微分方程.判断为恰当微分方程后，则可用如下解法:二 y 二 x设u(x, y) 是原方程的解，则 =m ,所以设u = fm(x, y)dx +v(y), ；x则卫=n二 yajm(x，y)dx+v(y) 一(fm(x，y)dx)-7；:y.誓，所以F( m(x, y)dx) dv(y) n 二:ydy,由此

6、v(y)= fnIm(x, y)dxdy ,由此可解得rr二 yu = Jm(x, y)dx + J nJm(x, y)dxdy ,所以原方程的通解为：:ym(x, y)dx - i n - m(x, y)dxdy = c, c 为任意常数。：:y首先可以求得其精确解为：y=exp(sin(x+2)> > x=-2:0.1:2;> > y= exp(sin(x+2)> > plot(x,y,'r.-');>> Data=x',y'3Columns 1 through 41.000000000000001.34382

7、524373165Columns 5 through 81.476121946445731.90449653438673Columns 9 through 122.049008650164272.43807150515633Columns 13 through 162.539682532380782.71148101768216Columns 17 through 202.717123008431282.57616043684702Columns 21 through 241.104986830331691.615146296442082.188741912604612.62100592628

8、6702.695718599203821.219778556000621.758818845766992.319776824715852.679016447572712.64811384739078-1.90000000000000 1.104986830331692.10792744704554Columns 25 through 281.964942888556771.819336991081061.674477827371160.90000000000000 1.27029521881156#1.53323499677732Columns 29 through 321.397923819

9、945001.270295218811561.151562836514531.04245724559826Columns 33 through 360.943296953050410.854066948581020.774497302497390.70413637458179Columns 37 through 400.642415207750370.588701425856340.542342319593710.50269776253737Column 410.46916418587400Data =-2.00000000000000 1.00000000000000-1.800000000

10、00000 1.21977855600062-1.70000000000000 1.34382524373165-1.60000000000000 1.47612194644573-1.50000000000000 1.61514629644208-1.40000000000000 1.75881884576699-1.30000000000000 1.90449653438673-1.20000000000000 2.04900865016427-1.10000000000000 2.18874191260461-1.00000000000000 2.31977682471585-0.900

11、00000000000 2.43807150515633-0.80000000000000 2.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.10000000000000 2.57616043684702

12、0 2.482577728015000.10000000000000 2.370757126170310.20000000000000 2.244530580577550.30000000000000 2.107927447045540.40000000000000 1.964942888556770.50000000000000 1.819336991081060.60000000000000 1.674477827371160.70000000000000 1.533234996777320.80000000000000 1.39792381994500有唯一解。则由欧拉法求初值问题(1)

13、, (2)的数值解的差分方程为:1.100000000000001.200000000000001.300000000000001.400000000000001.500000000000001.600000000000001.700000000000001.800000000000001.900000000000002.000000000000001.042457245598260.943296953050410.854066948581020.774497302497390.704136374581790.642415207750370.588701425856340.5423423195

14、93710.502697762537370.46916418587400#-2-1.5-1-0.500.511.525.1.2用欧拉法求解设常微分方程的初始问题dy = f(x, y),a 一 x b,(1) dxy(a) = .(2)yn 1 = Vn hf (Xn,yn), yi =.程序如下：建立函数文fl.mfunction x,y=f1(fun,x_span,y0,h)x=x_span(1):h:x_span(2);y(i)=y0;for n=1:length(x)-1y(n+1)=y(n)+h*feval(fun,x(n),y(n);endx=x'y=y'在MATL

15、AB输入以下程序：> > clear all> > fun=inline(' y*cos(x+2)');> > x,y1=f1(fun,-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.40000000000000-1.3000000000000

16、01.000000000000001.100000000000001.209450458180581.327984655342341.454851875167081.588852606593921.728287540690011.87092926670362-1.20000000000000 2.014025829963631.100000000000001.000000000000000.900000000000000.800000000000000.700000000000000.600000000000000.500000000000000.400000000000000.3000000

17、00000000.200000000000000.1000000000000000.100000000000000.200000000000000.300000000000000.400000000000000.500000000000000.600000000000000.700000000000000.800000000000000.900000000000001.000000000000001.100000000000001.200000000000001.300000000000001.400000000000001.500000000000001.600000000000002.15

18、4344360817052.288260553794212.411895799158422.521298457136512.612659661865862.682548001780242.728142503735772.747440620382272.739418225015642.704122329429032.642684103673772.557248883750392.450829780426752.327100593658172.190150463724832.044225990027361.893486050208131.741790624182991.59253854452440

19、1.448561571205111.312074863782761.184677883555021.067395661994220.960748409479460.864837397675810.779436454229260.704080678711320.638146572714181.700000000000000.580920241720559-1.700000000000001.34289777810139111.90000000000000 0.489600406402422.00000000000000 0.4540587312867232.521.510.50-2-1.5-1-

20、0.500.511.525.1.3用改进欧拉法求解：计算公式为：y0）1 =yn hf （Xn,yn）,ynkj =yn +f（Xn,yn） + f （Xn+ynki）,V。二 ,n=0,1,2,.,k =0,1,2,.,即先用欧拉法得（Xn,yn）,进而由（3）的第一式得初始近似值ynM,然后再用（3） 的第二式进行迭代，反复改进这个近似值，直到|二）-yk书|君（名为所允许的 误差）为止，并把取作为y（Xn+）的近似值yn+,这个方法就称为改进欧拉法。 通常称（3）为预报校正公式，其中第一式称为预报公式，第二式称为校正公式。 这个公式还可以写为：（下文改进欧拉法计算就以下面为准）yn 1.

21、二yn 2(kl k2),kl =f(Xn，yn),k2 =f (Xn h, yn hki), yo 二 J程序如下：建立函数文件f2.mfunction x,y=f2(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);y(n+1)=y(n)+h*k1;k2=feval(fun,x(n+1),y(n+1);y(n+1)=y(n)+h*(k1+k2)/2;endx=x'y=y'在MATLAB输入以下程序：> > clear all>

22、> fun=inline(' y*cos(x+2)');> > x,y2=f2(fun,-2,2,1,0.1);> > x,y2> > plot(x,y2,'b+-')结果及其图象：ans =-2.00000000000000-1.90000000000000-1.800000000000001.000000000000001.104725229090291.219207229363991.400000000000000.774969822614771.500000000000001.400000000000001.30

23、0000000000001.200000000000001.100000000000001.000000000000000.900000000000000.800000000000000.700000000000000.600000000000000.500000000000000.400000000000000.300000000000000.200000000000000.1000000000000000.100000000000000.200000000000000.300000000000000.400000000000000.500000000000000.6000000000000

24、00.700000000000000.800000000000000.900000000000001.000000000000001.100000000000001.200000000000001.613388617478771.756604945660731.901814952758942.045861837217642.185146633762432.315763555724452.433682968959042.534971671734682.616033679010132.673849667845002.706190767907232.711783264044682.690405219

25、408782.642903530441302.571129346284082.477799560853782.366300571501322.240456332742542.104285124930581.961768315515181.816650280351961.672282601880951.531518885262901.396660163895361.269445745893071.151080940502511.042291474098520.943394336221441.500000000000001.600000000000001.700000000000001.80000

26、0000000001.900000000000002.000000000000000.704729719126620.643092734337550.589432936264430.543103926871230.503471430640500.46993706594350132.5+十4-七七Hr+1.50.5-2-1.5-1-0.5七4-七00.511.525.1.4用4阶龙格一库塔求解标准的四阶龙格-库塔公式(亦称为经典的四阶龙格-库塔公式)hyn 1- =yn(k1 2k22k3 kJ,62k1 =f (xn,yn),1 hk2 uf (Xn 二 h, yn -k1), 22,一1.h

27、,、k3f (Xn-h,yn-k2),2 2k4 = f (Xn h, yn hka), y0 = ,n =0,1,2,.-1.500000000000001.61514577980337程序如下：建立函数文件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=feval(fun,x(n)+h/2,y(n)+h/2*k2);k4=feval(fun,x

28、(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 =1.400000000000000.77449743625256#-2.00000000000000-1.90

29、000000000000-1.80000000000000-1.70000000000000-1.600000000000001.000000000000001.104986745696811.219778373518281.343824953625391.476121543347421.400000000000001.300000000000001.200000000000001.100000000000001.000000000000000.900000000000000.800000000000000.700000000000000.600000000000000.50000000000

30、0000.400000000000000.300000000000000.200000000000000.1000000000000000.100000000000000.200000000000000.300000000000000.400000000000000.500000000000000.600000000000000.700000000000000.800000000000000.900000000000001.000000000000001.100000000000001.200000000000001.300000000000001.758818219617671.904495

31、806501982.049007830855662.188741013447442.319775857524332.438070481408852.539681463116192.621004822310582.679015319711092.711479876846582.717121865404832.695717464250942.648112729936422.576159345482252.482576670951542.370756112042032.244529619279342.107926550209661.964942069343181.819336263174131.67

32、4477203345411.533234486212941.397923427764521.270294944247731.151562672942311.042457181239250.943296972359340.854067034002241.50000000000000 0.704136540203991.600000000000001.700000000000001.800000000000001.900000000000002.000000000000000.642415391183080.588701616051610.542342508643010.5026979454353

33、60.469164360045032.51.50.501-2-1.5rr-1-0.500.511.525.1.5问题讨论与分析由以上数值分析结果绘制表格:精确解xi yi-2.00 1.000000-1.90 1.104987欧拉方法yi 误差1.000000 0.0000001.100000 0.004987改进的欧拉方法yi 误差1.000000 0.0000001.104725 0.000262四阶龙格-库塔方法yi误差1.0000000.000000-1.80 1.219779-1.70 1.343825-1.60 1.476122-1.50 1.615146-1.40 1.75881

34、9-1.30 1.904497-1.20 2.049009-1.10 2.1887421.2094501.3279851.4548521.5888531.7282881.8709292.0140262.1543440.0103280.0158410.0212700.0262940.0305310.0335670.0349830.0343981.2192071.3428981.4747971.6133891.7566051.9018152.0458622.1851470.0005710.0009270.0013250.0017580.0022140.0026820.0031470.0035951

35、.1049871.2197781.3438251.4761221.6151461.7588181.9044962.0490082.1887410.0000000.0000000.0000000.0000000.0000010.0000010.0000010.0000010.00000119-1.00 2.319777-0.90 2.438072-0.80 2.539683-0.70 2.621006-0.60 2.679016-0.50 2.711481-0.40 2.717123-0.30 2.695719-0.20 2.648114-0.10 2.5761600.002.4825780.1

36、02.3707570.202.2445310.302.1079270.401.9649430.501.8193370.601.6744780.701.5332350.801.3979240.901.2702951.001.1515631.101.0424571.200.9432971.300.8540671.400.7744971.500.7041361.600.6424151.700.5887011.800.5423421.900.5026972.000.4691642.288261 0.0315162.411896 0.0261762.521298 0.0183842.612660 0.0

37、083462.682548 0.0035322.728143 0.0166612.747441 0.0303182.739418 0.0437002.704122 0.0560082.642684 0.0665242.557249 0.0746712.450830 0.0800732.327101 0.0825702.190150 0.0822232.044226 0.0792831.893486 0.0741491.741791 0.0673131.592539 0.0593041.448562 0.0506381.312075 0.0417801.184678 0.0331151.0673

38、96 0.0249380.960748 0.0174510.864837 0.0107700.779436 0.0049390.704081 0.0000560.638147 0.0042690.580920 0.0077810.531652 0.0106900.489600 0.0130970.454058 0.0151062.315764 0.0040132.433683 0.0043892.534972 0.0047112.616034 0.0049722.673850 0.0051672.706191 0.0052902.711783 0.0053402.690405 0.005313

39、2.642904 0.0052102.571129 0.0050312.477800 0.0047782.366301 0.0044572.240456 0.0040742.104285 0.0036421.961768 0.0031751.816650 0.0026871.672283 0.0021951.531519 0.0017161.396660 0.0012641.269446 0.0008491.151081 0.0004821.042291 0.0001660.943394 0.0000970.854376 0.0003090.774970 0.0004730.704730 0.

40、0005930.643093 0.0006780.589433 0.0007320.543104 0.0007620.503471 0.0007740.469937 0.0007732.319776 2.438070 2.539681 2.621005 2.679015 2.711480 2.717122 2.695717 2.648113 2.576159 2.482577 2.370756 2.244530 2.107927 1.964942 1.819336 1.674477 1.533234 1.397923 1.270295 1.151563 1.042457 0.943297 0.

41、854067 0.774497 0.704137 0.642415 0.588702 0.542343 0.502697 0.4691640.0000010.0000010.0000010.0000010.0000010.0000010.0000010.0000010.0000010.0000010.0000010.0000010.0000010.0000010.0000010.0000010.0000010.0000010.0000000.0000000.0000000.0000000.0000000.0000000.0000000.0000000.0000000.0000000.00000

42、00.0000000.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-')精确解与欧拉方程比较3精确解与4阶龙格一库塔比较2100.51.520 1L1-2-1.5-1-0.5精确解与改进后的欧拉方程比较3由

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

44、r n=1:length(x)-1y(n+1)=y(n)+h*feval(fun,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.

45、2145-1.75001.2740-1.70001.3357-1.65001.3995-1.60001.4653-1.55001.5327-1.50001.6018-1.45001.6720-1.40001.7433-1.35001.8153-1.30001.8875-1.25001.9597-1.20002.0314-1.15002.1021-1.10002.1715-1.05002.2390-1.00002.304123-0.90002.4252-0.85002.4803-0.80002.5309-0.75002.5768-0.70002.6174-0.65002.6524-0.60002

46、.6814-0.55002.7042-0.50002.7205-0.45002.7301-0.40002.7330-0.35002.7290-0.30002.7182-0.25002.7007-0.20002.6766-0.15002.6462-0.10002.6097-0.05002.567602.52002.82 62.42.221.01.S1.41.2 1-2-1.8/ 614-1.2-10 80,6£一42 D取步长h=0.02程序如下：建立函数文f5.mfunction x,y=f5(fun,x_span,y0,h)x=x_span(1):h:x_span(2);y(i)=

47、y0;for n=1:length(x)-1y(n+1)=y(n)+h*feval(fun,x(n),y(n);endx=x'y=y'在MATLAB输入以下程序：>> clear all>> fun=inline(' y*cos(x+2)');>> x,y5=f5(fun,-2,0,1,0.02);>> x,y5>> plot(x,y5,'r+-')25-2.00001.0000-1.98001.0200-1.96001.0404-1.94001.0612-1.92001.0824-1.

48、90001.1040-1.88001.1259-1.86001.1483-1.84001.1710-1.82001.1941-1.80001.2176-1.78001.2415-1.76001.2657-1.74001.2903-1.72001.3153-1.70001.3405-1.68001.3662-1.66001.3921-1.64001.4183-1.62001.4449-1.60001.4717-1.58001.4988-1.56001.5262-1.54001.5538-1.52001.5817-1.50001.6097-1.48001.6380-1.46001.6664-1.4

49、4001.6950#-1.42001.7237-1.40001.7526-1.38001.7815-1.36001.8105-1.34001.8395-1.32001.8686-1.30001.8977-1.28001.9267-1.26001.9556-1.24001.9845-1.22002.0133-1.20002.0419-1.18002.0704-1.16002.0986-1.14002.1266-1.12002.1544-1.10002.1818-1.08002.2090-1.06002.2357-1.04002.2621-1.02002.2881-1.00002.3135-0.9

50、8002.3385-0.96002.3630-0.94002.3869-0.92002.4103-0.90002.4330-0.88002.4551-0.86002.4765-0.82002.5171-0.80002.5363-0.78002.5547-0.76002.5722-0.74002.5889-0.72002.6048-0.70002.6197-0.68002.6337-0.66002.6468-0.64002.6589-0.62002.6700-0.60002.6801-0.58002.6893-0.56002.6973-0.54002.7044-0.52002.7104-0.50

51、002.7153-0.48002.7191-0.46002.7219-0.44002.7235-0.42002.7241-0.40002.7236-0.38002.7220-0.36002.7194-0.34002.7156-0.32002.7108-0.30002.7049-0.28002.6979-0.26002.6899-0.24002.6808#-0.22002.6707-0.20002.6596-0.18002.6475-0.16002.6345-0.14002.6205-0.12002.6055-0.10002.5897-0.08002.5729-0.06002.5553-0.04002.5369-0.0