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

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

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

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

4、社会上与众不同的人的比例， xi (t)表示在时刻 t人口 中与众不同的人的数量。1)假定 p(0) 0.01,b 0.02和 r 0.1，当步长为 h 1年时，求从 t 0到t 50解 p(t)的近 似值，并作出近似解的曲线图形。2)精确求出微分方程的解 p(t) ，并将你当 t 50 时在分题 (b)中得到的结果与此时的精确值 进行比较。MATLAB 相关函数 】求微分方程的解析解及其数值的代入dsolve(egn1， egn2 ， L x ) subs (expr, x,y, x1,y1, )其中egn i 表示第 i 个方程， x 表示微分方程中的自变量，默认时自变量为t 。subs

5、命令中的 expr、 x、 y 为符合型表达式， x、y 分别用数值 x1、x2 代入。 syms x y z subs(x+y+z,x,y,z,1,2,3) ans =6 syms x subs(x2,x,2) ans =4s=dsolve(Dy 1 y 2, y(0) 1, x )ans =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

6、 (var1, var2, ) 。 f=inline(x+3*y,x,y)Inline function: f(x,y) = x+3*y f(2,3) ans =114， 5 阶龙格库塔方法求解微分方程数值解t,x=ode45(f,ts,x0,options)其中 f 是由待解方程写成的 m 文件名； x0 为函数的初值； t,x 分别为输出的自变量和函 数值(列向量 )，t 的步长是程序根据误差限自动选定的。若ts=t0,t1,t2, ,tf ，则输出在自变量指定值，等步长时用 ts=t0:k:tf ，输出在等分点； options 用于设定误差限 (可以缺省，缺省时设定为相对误差 10 3

7、 ，绝对误差 10 6 )，程序为：options=odeset( reltol ,rt, abstol,at)，这里 rt,at 分别为设定的相对误差和绝对误差。常用 选项见下表。选项名功能可选值省缺值AbsTol设定绝对误差正数1e 6RelTol设定相对误差正数1e 3InitialStep设定初始步长正数自动MaxStep设定步长上界正数tspan 10MaxOrder设定 ode15s 的最高阶数1,2,3,4,55Stats显示计算成本统计on,offoffBDF设定 ode15s 是否用反向差分on,offoff例： 解微分方程y y 2 t 0 t 4 yy(0) 1在命令窗口

8、执行odefun = inline (y 2* t y , t , y ) ；t,y ode45(odefun,0,4,1) ；t,yans =01.00000.05021.04900.10051.09590.15071.14083.85072.95033.90052.96723.95022.98394.00003.0006plot( t , y ,o-,)ode45(odefun,0,4,1)%解函数图形表示%不用输出变量，则直接输出图形t,y ode45(odefun,0 : 4,1) ； t,yans =01.00001.00001.73212.00002.23613.00002.645

9、84.00003.0006. 操作方法与实验步骤(包括实验数据记录和处理)2-1 编程编写用向前欧拉公式和改进欧拉公式求微分方程数值解的 Matlab 程序，问题如下： 在区间 a,b 内 (N 1)个等距点处，逼近下列初值问题的解，并对程序的每一句添上注释语 句。y f(x,y)a x b y(a) y0Euler 法 y=euler(a,b,n,y0,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;

10、for i=1:n;y(i+1)=y(i)+h*f(x(i),y(i);endplot(x,y)hold on% 求微分方程的精确解 x1=linspace(a,b,100);精确解为 s=dsolve(f1,b1,x) syms x y1=zeros(1,100);fori=1:100y1(i)=subs(s,x,x1(i);endplot(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;fori=1:n;T1

11、=f(x(i),y(i);T2=f(x(i+1),y(i)+h*T1);y(i+1)=y(i)+(h/ 2)*(T1+T2);endplot(x,y) hold on% 求微分方程的精确解 x1=linspace(a,b,100);精确解为 s=dsolve(f1,b1,x) syms x y1=zeros(1,100);for i=1:100 y1(i)=subs(s,x,x1(i);endplot(x1,y1,r)title( 红色代表精确解 )2-2 分析应用题假设等分区间数 n 100 ，用欧拉法和改进欧拉法在区间 ty (t) y(t) 20 y(0) 10并作出解的曲线图形， 同时

12、将方程的解析解也画在同一张图上， 的精度。(1)向前欧拉法 euler(0,10,100,10,inline(y-20,x,y),Dy=y-20,y(0)=10)ans =精确解为s =20 - 10*exp(x)ans =1.0e+005 *Columns 1 through 80.0001 0.0001 0.0001 0.0001 0.00010,10 内求解初值问题并作比较， 分析这两种方法0.0000 0.00000.0000Columns 9 through 16-0.0000 -0.0000-0.0001-0.0001-0.0001-0.0001-0.0002-0.0002Colu

13、mns 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.0012-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.007

14、8-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.0227-0.0250-0.0275-0.0302-0.0333-0.0366-0.0403Columns 65 through 72-0.0444 -0.0488-0.0537-0.0591-0.0651-0.0716-0.0788-0.0867Columns 73 through 80-0.0954 -0.1049 -0.1154 -0.1270 -0.13

15、97 -0.1537 -0.1691-0.1860Columns 81 through 88-0.2046 -0.2251-0.2477-0.2724-0.2997-0.3297-0.3627-0.3990Columns 89 through 96-0.4389 -0.4828-0.5311-0.5842-0.6427-0.7070-0.7777-0.8555Columns 97 through 101-0.9410 -1.0352-1.1387-1.2526-1.37792）改进欧拉法 eulerpro(0,10,100,10,inline(y-20,x,y),Dy=y-20,y(0)=10

16、) ans =精确解为20 - 10*exp(x) ans =1.0e+005 *Columns 1 through 80.0001 0.0001 0.00010.00010.00010.00000.0000-0.0000Columns 9 through 16-0.0000 -0.0000 -0.0001-0.0001-0.0001-0.0002-0.0002-0.0002Columns 17 through 24-0.0003 -0.0003 -0.0004-0.0005-0.0005-0.0006-0.0007-0.0008Columns 25 through 32-0.0009 -0.

17、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-0.0042-0.0047Columns 41 through 48-0.0052 -0.0058 -0.0064-0.0071-0.0079-0.0087-0.0097-0.0107Columns 49 through 56-0.0119 -0.0131 -0.0145-0.0161-0.0178-0.0197-0.0218-0.0241Columns 57 throu

18、gh 64-0.0266 -0.0294 -0.0325-0.0360-0.0398-0.0440-0.0486-0.0537Columns 65 through 72-0.0594 -0.0656 -0.0726-0.0802-0.0886-0.0980-0.1083-0.1197Columns 73 through 80-0.1323 -0.1462 -0.1615-0.1785-0.1973-0.2180-0.2409-0.2663Columns 81 through 88-0.2942 -0.3251 -0.3593-0.3971-0.4388-0.4849-0.5358-0.5921

19、Columns 89 through 96-0.6543 -0.7230 -0.7989-0.8828-0.9755-1.0780 -1.1912-1.3163Columns 97 through 101-1.4545 -1.6073 -1.7760 -1.9626 -2.1686改进欧拉法的精度比向前欧拉法更高。2-3 分析应用题 用以下三种不同的方法求下述微分方程的数值解，取 h 10y y 2x 0 x 1y(0) 1画出解的图形，与精确值比较并进行分析。1)欧拉法；2)改进欧拉法；2-4 分析应用题 考虑一个涉及到社会上与众不同的人的繁衍问题模型。假设在时刻t (单位为年 )，社会上有

20、人口 x(t) 人，又假设所有与众不同的人与别的与众不同的人结婚后所生后代也是与众 不同的人。 而固定比例为 r 的所有其他的后代也是与众不同的人。 如果对所有人来说出生率 假定为常数 b ，又如果普通的人和与众不同的人的婚配是任意的， 则此问题可以用微分方程 表示为：dp(t)rb(1 p(t) dt其中变量 p(t) xi(t) x(t)表示在时刻 t 社会上与众不同的人的比例， xi (t)表示在时刻 t人口 中与众不同的人的数量。1)假定 p(0) 0.01,b 0.02和 r 0.1，当步长为 h 1年时，求从 t 0到t 50解 p(t)的近 似值，并作出近似解的曲线图形。2)精确求出微分方程的解 p(t)，并将你当 t 50 时在分题 (b)中得到的结果与此时的精确值 进行比较。1) euler(0,50,50,0.01,inline(0.