实验报告七常微分方程初值问题的数值解法

1、浙江大学城市学院实验报告课程名称数值计算方法实验项目名称常微分方程初值问题的数值解法实验成绩 指导老师(签名) 日期 2015/12/16一、实验目的与要求1. 用Matlab软件掌握求微分方程数值解的欧拉方法与龙格-库塔方法；2. 通过实例学习用微分方程模型解决简化的实际问题。二、实验内容与原理编程题2-1要求写出Matlab源程序(m文件),并有适当的注释语句;分析应用题2-2,2-3,2-4,2-5 要求将问题的分析过程、Matlab源程序与运行结果与结果的解释、算法的分析写在实验报 告上。2-1编程编写用向前欧拉公式与改进欧拉公式求微分方程数值解的Matlab程序，问题如下：在区间a,

2、b内(N 1)个等距点处，逼近下列初值问题的解，并对程序的每一句添上注释语句。y f (x, y) a x b y(a) yEuler 法 y=euler(a,b,n,y0,f,f1,b1)改进 Euler 法 y=eulerpro(a,b,n,y0,f,f1,b1)2-2分析应用题假设等分区间数 n 100,用欧拉法与改进欧拉法在区间t 0,10内求解初值问题y(t) y(t) 20y(0)10并作出解的曲线图形，同时将方程的解析解也画在同一张图上，并作比较，分析这两种方法的精度。2-3分析应用题用以下三种不同的方法求下述微分方程的数值解，取h 10y y 2x 0 x 1y(0) 1画出解

3、的图形，与精确值比较并进行分析。1) 欧拉法；2) 改进欧拉法；3) 龙格库塔方法；2-4分析应用题考虑一个涉及到社会上与众不同的人的繁衍问题模型。假设在时刻t(单位为年),社会上有人口 x(t)人,又假设所有与众不同的人与别的与众不同的人结婚后所生后代也就是与 众不同的人。而固定比例为r的所有其她的后代也就是与众不同的人。如果对所有人来说出生率假定为常数b ,又如果普通的人与与众不同的人的婚配就是任意的，则此问题可以用微分方程表示为：dP(t)rb(1 p(t)dt其中变量p(t)Xi(t).：x(t)表示在时刻t社会上与众不同的人的比例 ，Xi(t)表示在时刻t人口中与众不同的人的数量。1

4、) 假定p(0)0.01,b0.02与r 0.1,当步长为h 1年时，求从t 0到t 50解p(t)的近似值,并作出近似解的曲线图形。2) 精确求出微分方程的解p(t),并将您当t 50时在分题(b)中得到的结果与此时的精确值进行比较。【MATLAB相关函数】求微分方程的解析解及其数值的代入dsolve( egn1 dgn2,LX )subs (expr, x,y, , x1,y1, )其中egni 表示第i个方程，x 表示微分方程中的自变量，默认时自变量为t。subs命令中的expr、x、y为符合型表达式，x、y分别用数值x1、x2代入。 syms x y z subs(x+y+z,x,y,

5、z,1,2,3)ans =6 syms x subs(xA2,x,2)ans =4s=dsolve( Dy 1 y 2 y(0)1, xans =tan(x 1 4 pi) syms x subs(s,x,2)ans =-0、 3721右端函数f (x, y)的自动生成f= inline( expr var1 Var2 )其中expr表示函数的表达式，var1, Var2 表示函数表达式中的变量，运行该函数 住成一个新的函数表达式为f (var1, var2,)。 f=i nlin e(x+3*y,x,y)f =Inline fun ctio n:f(x,y) = x+3*y f(2,3)an

6、s =114,5阶龙格-库塔方法求解微分方程数值解t,x=ode45(f,ts,x0,optio ns)其中f就是由待解方程写成的m文件名；x0为函数的初值;t,x分别为输出的自变量与函数值(列向量),t的步长就是程序根据误差限自动选定的。若ts=t0,t1,t2,tf,则输出在自变量指定值，等步长时用ts=tO:k:tf,输出在等分点;options用于设定误差限(可以缺省，缺省时设定为相对误差10 3,绝对误差10 6),程序差与绝对误差。常用选项见下表。选项名功能可选值省缺值AbsTol设定绝对误差正数1e 6RelTol设定相对误差正数1e 3In itialStep设疋初始步长正数自

7、动MaxStep2 rc 宀 rH x . ffl 设定步长上界正数tspa n/10MaxOrder设定ode15s的最高阶数123,4,55Stats显示计算成本统计on ,offoff BDF设定ode15s就是否用反向差 分on ,offoff例:解微分方程y y 20 t 4yy(0)1在命令窗口执行odefun = inline ( 2* t y , t, )；t, yode45(odefun,0,4,1);t,yans =plot( t , y, 6-)ode45(odefu n,0,4,1)t, y ode45(odefun,0 : 4,1)；ans =%解函数图形表示%不用输

8、出变量，则直接输出图形t,y0、005021、 00001、04900、10051、09590、15071、14083、85072、95033、90052、96723、95022、98394、00003、000601、 00001、00001、73212、00002、23613、00002、64584、00003、0006操作方法与实验步骤(包括实验数据记录与处理)2-1编程编写用向前欧拉公式与改进欧拉公式求微分方程数值解的Matlab程序，问题如下：在区间a,b内(N 1)个等距点处，逼近下列初值问题的解，并对程序的每一句添上注释语 句。y f (x, y) a x b y(a)yEuler

9、 法y=euler(a,b,n,yO,f,f1,b1)改进 Euler 法 y=eulerpro(a,b,n,y0,f,f1,b1)Euler 法y=euler(a,b,n,y0,f,f1,b1) y=zeros(1,n+1); y(1)=y0;h=(b-a)/n; x=a:h:b; for i=1:n;y(i+1)=y(i)+h*f(x(i),y(i); end plot(x,y) hold on% 求微分方程的精确解 x1=linspace(a,b,100); 精确解为 s=dsolve(f1,b1,x) syms x y1=zeros(1,100); fori=1:100y1(i)=su

10、bs(s,x,x1(i); end plot(x1,y1,r) title( 红色代表精确解 ) 改进 Euler 法 y=eulerpro(a,b,n,y0,f,f1,b1) % 求微分方程的数值解 y=zeros(1,n+1);y(1)=y0; h=(b-a)/n; x=a:h:b; for i=1:n; T1=f(x(i),y(i); T2=f(x(i+1),y(i)+h*T1); y(i+1)=y(i)+(h/2)*(T1+T2); end plot(x,y) hold on% 求微分方程的精确解 x1=linspace(a,b,100);精确解为 s=dsolve(f1,b1,x)

11、syms x y1=zeros(1,100); for i=1:100y1(i)=subs(s,x,x1(i);endplot(x1,y1,r)title( 红色代表精确解 )2-2 分析应用题 假设等分区间数 n 100, 用欧拉法与改进欧拉法在区间 t 0,10 内求解初值问题y(t) y(t) 20 y(0) 10并作出解的曲线图形 , 同时将方程的解析解也画在同一张图上 , 并作比较 , 分析这两种方法的 精度。(1)向前欧拉法 euler(0,10,100,10,inline(y-20,x,y),Dy=y-20,y(0)=10)ans =精确解为s =20 - 10*exp(x) a

12、ns =1、0e+005 *Columns 1 through 80、 00010、00010、00010、00010、00010、00000、00000、0000Columns 9 through 16-0、 0000 -0、 0000-0、0001-0、0001-0、0001-0、0001-0、 0002-0、0002Columns 17 through 24-0、 0003 -0、 0003-0、0004-0、0004-0、0005-0、0005-0、 0006-0、0007Columns 25 through 32-0、 0008 -0、 0009-0、0010-0、0011-0、00

13、12-0、0014-0、 0015-0、0017Columns 33 through 40-0、 0019 -0、 0021-0、0024-0、0026-0、0029-0、0032-0、 0035-0、0039Columns 41 through 48-0、 0043 -0、 0048-0、0053-0、0058-0、0064-0、0071-0、 0078-0、0086Columns 49 through 56-0、 0095 -0、 0105-0、0115-0、0127-0、0140-0、0154-0、 0170-0、0187Columns 57 through 64-0、 0206 -0、

14、 0227-0、0250-0、0275-0、0302-0、0333-0、 0366-0、0403Colu mns 65 through 72-0、-0、 0444-0、 0488-0、0537-0、0867Columns 73 through 80-0、 0954-0、 1049-0、1154-0、-0、1860Columns 81 through 88-0、 2046-0、 2251-0、2477-0、-0、3990Columns 89 through 96-0、 4389-0、 4828-0、5311-0、-0、8555Columns 97 through 101-0、9410-1、035

15、2-1、1387-1、2526-1、37790591-0、0651-0、0716-0、07881270-0、1397-0、1537-0、16912724-0、2997-0、3297-0、36275842-0、6427-0、7070-0、7777改进欧拉法g *- eulerpro(0,10,100,10,i nlin e(y-20,x,y Hattd吐.(2)ans =20 - 10*exp(x)精确解为ans =1、 0e+005*Columns 1 through 8b Q0000 0、)0000、0001 0、0001 0、0、 0001 0、 00010、0001-0、Columns

16、9 through 16-0、0000 -15-0、 0000-0、0001-0、 0001-0、 0001-0、 0002-0、 00020002Colu mns 17 through 24-0、 00030、 0003-0、0004-0、 0005-0、 0005-0、 0006-0、 00070000-0、-0、Colu mns 25gh 320008-0、0009-0、0010-0、0011-0、0013-0、0014-0、0016-0、0018-0、 0020Columns 33 through 40-0、0022-0、0025-0、0028-0、0031-0、0034-0、0038

17、-0、0042-0、0047Colu mns 41 through 48-0、 0052-0、 0058-0、0064-0、0071-0、0079-0、0087-0、0097-0、0107Colu mns 49 through 56-0、0119-0、0131-0、0145-0、0161-0、0178-0、0197-0、0218-0、Columns 57 through 64-0、 0266-0、 0294-0、 0325-0、0537Columns 65 through 72-0、 0594-0、 0656-0、 0726-0、1197Columns 73 through 80-0、1323

18、-0、1462-0、1615-0、2663Columns 81 through 88-0、 2942-0、 3251-0、 3593-0、5921Columns 89 through 96-0、 6543-0、 7230-0、 7989-0、3163Colu mns 97 through 101-1、4545-1、6073-1、7760-1、9626-2、 16860360-0、0398-0、0440-0、0486-0、0802-0、0886-0、0980-0、1083-0、1785-0、1973-0、2180-0、2409-0、3971-0、4388-0、4849-0、5358-0、8828

19、-0、9755-1、0780-1、1912-1、.ni2d改进欧拉法的精度比向前欧拉法更高。2-3分析应用题用以下三种不同的方法求下述微分方程的数值y y 2x 0 y(0) 1画出解的图形，与精确值比较并进行分析。1) 欧拉法；2) 改进欧拉法；2-4分析应用题U* tutBnairt luili心考虑一个涉及到社会上与众不同的人的繁衍问题模型。假设在时刻t(单位为年),社会上有人口 x(t)人,又假设所有与众不同的人与别的与众不同的人结婚后所生后代也就是与众不同的人。而固定比例为r的所有其她的后代也就是与众不同的人。如果对所有人来说出生率假定为常数 b ,又如果普通的人与与众不同的人的婚配

20、就是任意的，则此问题可以用微分方程表示为：dPrb(1 p(t)dt其中变量p(t) Xi(t)；X(t)表示在时刻t社会上与众不同的人的比例，Xi(t)表示在时刻t人口中与众不同的人的数量。1) 假定p(0)0.01,b0.02与r 0.1,当步长为h 1年时，求从t 0到t 50解p(t)的近似值,并作出近似解的曲线图形。2) 精确求出微分方程的解p(t),并将您当t 50时在分题(b)中得到的结果与此时的精确值进行比较。euler(0,50,50,0、01,inline(O、002-0、002*p,t,p),Dp=0、002-0、002*p,p(0)=0、0 1)ans =精确解为1 - 99/(10