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

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

2、区间内个等距点处，逼近下列初值问题的解，并对程序的每一句添上注释语句。 Euler法 y=euler(a,b,n,y0,f,f1,b1)改进Euler法 y=eulerpro(a,b,n,y0,f,f1,b1)2-2 分析应用题假设等分区间数，用欧拉法和改进欧拉法在区间内求解初值问题并作出解的曲线图形，同时将方程的解析解也画在同一张图上，并作比较，分析这两种方法的精度。2-3 分析应用题用以下三种不同的方法求下述微分方程的数值解，取 画出解的图形，与精确值比较并进行分析。 1）欧拉法； 2）改进欧拉法；3）龙格库塔方法； 2-4 分析应用题考虑一个涉及到社会上与众不同的人的繁衍问题模型。假设在

3、时刻(单位为年)，社会上有人口人，又假设所有与众不同的人与别的与众不同的人结婚后所生后代也是与众不同的人。而固定比例为的所有其他的后代也是与众不同的人。如果对所有人来说出生率假定为常数，又如果普通的人和与众不同的人的婚配是任意的，则此问题可以用微分方程表示为： 其中变量表示在时刻社会上与众不同的人的比例，表示在时刻人口中与众不同的人的数量。 1）假定和，当步长为年时，求从到解的近似值，并作出近似解的曲线图形。2）精确求出微分方程的解，并将你当时在分题(b)中得到的结果与此时的精确值进行比较。【MATLAB相关函数】n 求微分方程的解析解及其数值的代入 dsolve(egn1， egn2， )

4、subs (expr, x,y, x1,y1, )其中egn表示第个方程，表示微分方程中的自变量，默认时自变量为。subs命令中的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 = 4 s=dsolve(, , )ans = >> syms x>> subs(s,x,2)ans = -0.3721n 右端函数的自动生成 f

5、= inline(expr, var1, var2,) 其中expr表示函数的表达式，var1, var2 表示函数表达式中的变量，运行该函数，生成一个新的函数表达式为f (var1, var2, )。>> f=inline('x+3*y','x','y')f = Inline function: f(x,y) = x+3*y>> f(2,3)ans = 11n 4，5阶龙格库塔方法求解微分方程数值解t,x=ode45(f,ts,x0,options) 其中f是由待解方程写成的m文件名；x0为函数的初值；t,x分别为输出的

6、自变量和函数值(列向量)，t的步长是程序根据误差限自动选定的。若ts=t0,t1,t2,tf，则输出在自变量指定值，等步长时用ts=t0:k:tf，输出在等分点；options用于设定误差限(可以缺省，缺省时设定为相对误差，绝对误差)，程序为：options=odeset(reltol,rt,abstol,at)，这里rt,at分别为设定的相对误差和绝对误差。常用选项见下表。选项名功能可选值省缺值AbsTol设定绝对误差正数RelTol设定相对误差正数InitialStep设定初始步长正数自动MaxStep设定步长上界正数MaxOrder设定ode15s的最高阶数1,2,3,4,55Stats

7、显示计算成本统计on,offoffBDF设定ode15s是否用反向差分on,offoff例：解微分方程 在命令窗口执行 = (, , )； ； ans = 0 1.0000 0.0502 1.0490 0.1005 1.0959 0.1507 1.1408 3.8507 2.9503 3.9005 2.9672 3.9502 2.9839 4.0000 3.0006 plot(,o-,) %解函数图形表示 %不用输出变量，则直接输出图形 ； ans = 0 1.0000 1.0000 1.7321 2.0000 2.2361 3.0000 2.6458 4.0000 3.0006三. 操作方法

8、与实验步骤（包括实验数据记录和处理）2-1 编程编写用向前欧拉公式和改进欧拉公式求微分方程数值解的Matlab程序，问题如下：在区间内个等距点处，逼近下列初值问题的解，并对程序的每一句添上注释语句。 Euler法 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; for i=1:n; y(i+1)=y(i)+h*f(x(i),y(i); end plot(x,y) h

9、old on % 求微分方程的精确解 x1=linspace(a,b,100); '精确解为' s=dsolve(f1,b1,'x') syms xy1=zeros(1,100); for i=1:100 y1(i)=subs(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; T

10、1=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') syms x y1=zeros(1,100); for i=1:100 y1(i)=subs(s,x,x1(i); end plot(x1,y1,'r')title('红色代表精确解')2-2分析应用题假设等分区间数，用欧拉法和改进欧拉法在区

11、间内求解初值问题并作出解的曲线图形，同时将方程的解析解也画在同一张图上，并作比较，分析这两种方法的精度。(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 8 0.0001 0.0001 0.0001 0.0001 0.0001 0.0000 0.0000 0.0000 Column

12、s 9 through 16 -0.0000 -0.0000 -0.0001 -0.0001 -0.0001 -0.0001 -0.0002 -0.0002 Columns 17 through 24 -0.0003 -0.0003 -0.0004 -0.0004 -0.0005 -0.0005 -0.0006 -0.0007 Columns 25 through 32 -0.0008 -0.0009 -0.0010 -0.0011 -0.0012 -0.0014 -0.0015 -0.0017 Columns 33 through 40 -0.0019 -0.0021 -0.0024 -0.

13、0026 -0.0029 -0.0032 -0.0035 -0.0039 Columns 41 through 48 -0.0043 -0.0048 -0.0053 -0.0058 -0.0064 -0.0071 -0.0078 -0.0086 Columns 49 through 56 -0.0095 -0.0105 -0.0115 -0.0127 -0.0140 -0.0154 -0.0170 -0.0187 Columns 57 through 64 -0.0206 -0.0227 -0.0250 -0.0275 -0.0302 -0.0333 -0.0366 -0.0403 Colum

14、ns 65 through 72 -0.0444 -0.0488 -0.0537 -0.0591 -0.0651 -0.0716 -0.0788 -0.0867 Columns 73 through 80 -0.0954 -0.1049 -0.1154 -0.1270 -0.1397 -0.1537 -0.1691 -0.1860 Columns 81 through 88 -0.2046 -0.2251 -0.2477 -0.2724 -0.2997 -0.3297 -0.3627 -0.3990 Columns 89 through 96 -0.4389 -0.4828 -0.5311 -

15、0.5842 -0.6427 -0.7070 -0.7777 -0.8555 Columns 97 through 101 -0.9410 -1.0352 -1.1387 -1.2526 -1.3779（2） 改进欧拉法 >> eulerpro(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 8 0.0001 0

16、.0001 0.0001 0.0001 0.0001 0.0000 0.0000 -0.0000 Columns 9 through 16 -0.0000 -0.0000 -0.0001 -0.0001 -0.0001 -0.0002 -0.0002 -0.0002 Columns 17 through 24 -0.0003 -0.0003 -0.0004 -0.0005 -0.0005 -0.0006 -0.0007 -0.0008 Columns 25 through 32 -0.0009 -0.0010 -0.0011 -0.0013 -0.0014 -0.0016 -0.0018 -0

17、.0020 Columns 33 through 40 -0.0022 -0.0025 -0.0028 -0.0031 -0.0034 -0.0038 -0.0042 -0.0047 Columns 41 through 48 -0.0052 -0.0058 -0.0064 -0.0071 -0.0079 -0.0087 -0.0097 -0.0107 Columns 49 through 56 -0.0119 -0.0131 -0.0145 -0.0161 -0.0178 -0.0197 -0.0218 -0.0241 Columns 57 through 64 -0.0266 -0.029

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

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

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