实验报告七常微分方程初值问答的数值解法Word文件下载.docx

1、用以下三种不同的方法求下述微分方程的数值解，取画出解的图形，与精确值比较并进行分析。1）欧拉法；2）改进欧拉法；3）龙格库塔方法；2-4 分析应用题考虑一个涉及到社会上与众不同的人的繁衍问题模型。假设在时刻（单位为年），社会上有人口人，又假设所有与众不同的人与别的与众不同的人结婚后所生后代也是与众不同的人。而固定比例为的所有其他的后代也是与众不同的人。如果对所有人来说出生率假定为常数，又如果普通的人和与众不同的人的婚配是任意的，则此问题可以用微分方程表示为：其中变量表示在时刻社会上与众不同的人的比例，表示在时刻人口中与众不同的人的数量。1）假定和，当步长为年时，求从到解的近似值，并作出近似解的

2、曲线图形。2）精确求出微分方程的解，并将你当时在分题（b）中得到的结果与此时的精确值进行比较。【MATLAB相关函数】 求微分方程的解析解及其数值的代入 dsolve（egn1， egn2， ） 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 xx2,x,2） 4 s=dsolve（, , ）ans = subs（s,x,2） -0.3721 右端函数的自

3、动生成 f= inline（expr, var1, var2,） 其中expr表示函数的表达式，var1, var2 表示函数表达式中的变量，运行该函数，生成一个新的函数表达式为f （var1, var2, ）。 f=inline（x+3*y,xy）f = Inline function: f（x,y） = x+3*y f（2,3） 11 4，5阶龙格库塔方法求解微分方程数值解t,x=ode45（f,ts,x0,options） 其中f是由待解方程写成的m文件名；x0为函数的初值；t,x分别为输出的自变量和函数值（列向量），t的步长是程序根据误差限自动选定的。若ts=t0,t1,t2,tf，则

4、输出在自变量指定值，等步长时用ts=t0:k:tf，输出在等分点；options用于设定误差限（可以缺省，缺省时设定为相对误差，绝对误差），程序为：options=odeset（reltol,rt,abstol,at），这里rt,at分别为设定的相对误差和绝对误差。常用选项见下表。选项名功能可选值省缺值AbsTol设定绝对误差正数RelTol设定相对误差InitialStep设定初始步长自动MaxStep设定步长上界MaxOrder设定ode15s的最高阶数1,2,3,4,55Stats显示计算成本统计on,offoffBDF设定ode15s是否用反向差分例：解微分方程在命令窗口执行 = （,

5、 , ）； ； 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-,） %解函数图形表示 %不用输出变量，则直接输出图形 1.0000 1.7321 2.0000 2.2361 3.0000 2.6458三. 操作方法与实验步骤（包括实验数据记录和处理）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

6、 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,y1=zeros（1,100）;for i=1:100 y1（i）=subs（s,x,x1（i）;endplot（x1,y1,rtitle（红色代表精确解改进Euler法 y=eulerpro（a,b,n,y0,f,f1,b1） % 求微分方程的数值解 y=zeros（1,n+1）; y（1）=y0; h=（b-a）/n; x=a: for T1=f（x（i）,y（i）; T2=

7、f（x（i+1）,y（i）+h*T1）; y（i+1）=y（i）+（h/2）*（T1+T2）;% 求微分方程的精确解 y1=zeros（1,100）; for i=1:2-2分析应用题（1）向前欧拉法 euler（0,10,100,10,inline（y-20）,Dy=y-20y（0）=10） 精确解为s = 20 - 10*exp（x） 1.0e+005 * Columns 1 through 8 0.0001 0.0001 0.0001 0.0001 0.0001 0.0000 0.0000 0.0000 Columns 9 through 16 -0.0000 -0.0000 -0.00

8、01 -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.0026 -0.0029 -0.0032 -0.0035 -0.0039

9、 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 Columns 65 through 72 -0.0444 -0.0488 -0.

10、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 -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（ ans = 精确解为 s = 20 - 10*exp（x） 1.0e+00