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

-.z.浙江大学城市学院实验报告课程名称数值计算方法实验工程名称常微分方程初值问题的数值解法实验成绩指导教师〔签名〕日期2015/12/16一.实验目的和要求用Matlab软件掌握求微分方程数值解的欧拉方法和龙格－库塔方法；通过实例学习用微分方程模型解决简化的实际问题。二.实验内容和原理编程题2-1要求写出Matlab源程序(m文件)，并有适当的注释语句；分析应用题2-2，2-3，2-4，2-5要求将问题的分析过程、Matlab源程序和运行结果和结果的解释、算法的分析写在实验报告上。编程编写用向前欧拉公式和改良欧拉公式求微分方程数值解的Matlab程序，问题如下：在区间个等距点处，逼近以下初值问题的解，并对程序的每一句添上注释语句。Euler法y=euler(a,b,n,y0,f,f1,b1)改良Euler法y=eulerpro(a,b,n,y0,f,f1,b1)分析应用题假设等分区间数，用欧拉法和改良欧拉法在区间内求解初值问题并作出解的曲线图形，同时将方程的解析解也画在同一张图上，并作比拟，分析这两种方法的精度。分析应用题用以下三种不同的方法求下述微分方程的数值解，取画出解的图形，与准确值比拟并进展分析。1〕欧拉法；2〕改良欧拉法；3〕龙格－库塔方法；分析应用题考虑一个涉及到社会上与众不同的人的繁衍问题模型。假设在时刻(单位为年)，社会上有人口人，又假设所有与众不同的人与别的与众不同的人结婚后所生后代也是与众不同的人。而固定比例为的所有其他的后代也是与众不同的人。如果对所有人来说出生率假定为常数，又如果普通的人和与众不同的人的婚配是任意的，则此问题可以用微分方程表示为：其中变量表示在时刻社会上与众不同的人的比例，表示在时刻人口中与众不同的人的数量。1〕假定和，当步长为年时，求从到解的近似值，并作出近似解的曲线图形。2〕准确求出微分方程的解，并将你当时在分题(b)中得到的结果与此时的准确值进展比拟。【MATLAB相关函数】求微分方程的解析解及其数值的代入dsolve(‘egn1’，‘egn2’，‘’)subs(e*pr,{*,y,…},{*1,y1,…})其中‘egn’表示第个方程，‘’表示微分方程中的自变量，默认时自变量为。subs命令中的e*pr、*、y为符合型表达式，*、y分别用数值*1、*2代入。>>syms*yz>>subs('*+y+z',{*,y,z},{1,2,3})ans=6>>syms*>>subs('*^2',*,2)ans=4s=dsolve(‘’,‘’,‘’)ans=>>syms*>>subs(s,*,2)ans=-0.3721右端函数的自动生成f=inline(‘e*pr’,’var1’,‘var2’,其中’e*pr’表示函数的表达式，’var1’,‘var2’表示函数表达式中的变量，运行该函数，生成一个新的函数表达式为f(var1,var2,>>f=inline('*+3*y','*','y')f=Inlinefunction:f(*,y)=*+3*y>>f(2,3)ans=114，5阶龙格－库塔方法求解微分方程数值解[t,*]=ode45(f,ts,*0,options)其中f是由待解方程写成的m文件名；*0为函数的初值；t,*分别为输出的自变量和函数值(列向量)，t的步长是程序根据误差限自动选定的。假设ts=[t0,t1,t2,…,tf]，则输出在自变量指定值，等步长时用ts=t0:k:tf，输出在等分点；options用于设定误差限(可以缺省，缺省时设定为相对误差，绝对误差)，程序为：options=odeset(‘reltol’,rt,’abstol’,at)，这里rt,at分别为设定的相对误差和绝对误差。常用选项见下表。选项名功能可选值省缺值AbsTol设定绝对误差正数RelTol设定相对误差正数InitialStep设定初始步长正数自动Ma*Step设定步长上界正数Ma*Order设定ode15s的最高阶数1,2,3,4,55Stats显示计算本钱统计on,offoffBDF设定ode15s是否用反向差分on,offoff例：解微分方程在命令窗口执行=(‘’,‘’,‘’)；；ans=01.00000.05021.04900.10051.09590.15071.14083.85072.95033.90052.96723.95022.98394.00003.0006plot(,,‘o-’,)%解函数图形表示%不用输出变量，则直接输出图形；ans=01.00001.00001.73212.00002.23613.00002.64584.00003.0006三.操作方法与实验步骤〔包括实验数据记录和处理〕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;*=a:h:b;fori=1:n;y(i+1)=y(i)+h*f(*(i),y(i));endplot(*,y)holdon%求微分方程的准确解*1=linspace(a,b,100);'准确解为's=dsolve(f1,b1,'*')syms*y1=zeros(1,100);fori=1:100y1(i)=subs(s,*,*1(i));endplot(*1,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;*=a:h:b;fori=1:n;T1=f(*(i),y(i));T2=f(*(i+1),y(i)+h*T1);y(i+1)=y(i)+(h/2)*(T1+T2);endplot(*,y)holdon%求微分方程的准确解*1=linspace(a,b,100);'准确解为's=dsolve(f1,b1,'*')syms*y1=zeros(1,100);fori=1:100y1(i)=subs(s,*,*1(i));endplot(*1,y1,'r')title('红色代表准确解')2-2分析应用题假设等分区间数，用欧拉法和改良欧拉法在区间内求解初值问题并作出解的曲线图形，同时将方程的解析解也画在同一张图上，并作比拟，分析这两种方法的精度。(1)向前欧拉法>>euler(0,10,100,10,inline('y-20','*','y'),'Dy=y-20','y(0)=10')ans=准确解为s=20-10*e*p(*)ans=1.0e+005*Columns1through80.00010.00010.00010.00010.00010.00000.00000.0000Columns9through16-0.0000-0.0000-0.0001-0.0001-0.0001-0.0001-0.0002-0.0002Columns17through24-0.0003-0.0003-0.0004-0.0004-0.0005-0.0005-0.0006-0.0007Columns25through32-0.0008-0.0009-0.0010-0.0011-0.0012-0.0014-0.0015-0.0017Columns33through40-0.0019-0.0021-0.0024-0.0026-0.0029-0.0032-0.0035-0.0039Columns41through48-0.0043-0.0048-0.0053-0.0058-0.0064-0.0071-0.0078-0.0086Columns49through56-0.0095-0.0105-0.0115-0.0127-0.0140-0.0154-0.0170-0.0187Columns57through64-0.0206-0.0227-0.0250-0.0275-0.0302-0.0333-0.0366-0.0403Columns65through72-0.0444-0.0488-0.0537-0.0591-0.0651-0.0716-0.0788-0.0867Columns73through80-0.0954-0.1049-0.1154-0.1270-0.1397-0.1537-0.1691-0.1860Columns81through88-0.2046-0.2251-0.2477-0.2724-0.2997-0.3297-0.3627-0.3990Columns89through96-0.4389-0.4828-0.5311-0.5842-0.6427-0.7070-0.7777-0.8555Columns97through101-0.9410-1.0352-1.1387-1.2526-1.3779改良欧拉法>>eulerpro(0,10,100,10,inline('y-20','*','y'),'Dy=y-20','y(0)=10')ans=准确解为s=20-10*e*p(*)ans=1.0e+005*Columns1through80.00010.00010.00010.00010.00010.00000.0000-0.0000Columns9through16-0.0000-0.0000-0.0001-0.0001-0.0001-0.0002-0.0002-0.0002Columns17through24-0.0003-0.0003-0.0004-0.0005-0.0005-0.0006-0.0007-0.0008Columns25through32-0.0009-0.0010-0.0011-0.0013-0.0014-0.0016-0.0018-0.0020Columns33through40-0.0022-0.0025-0.0028-0.0031-0.0034-0.0038-0.0042-0.0047Columns41through48-0.0052-0.0058-0.0064-0.0071-0.0079-0.0087-0.0097-0.0107Columns49through56-0.0119-0.0131-0.0145-0.0161-0.0178-0.0197-0.0218-0.0241Columns57through64-0.0266-0.0294-0.0325-0.0360-0.0398-0.0440-0.0486-0.0537Columns65through72-0.0594-0.0656-0.0726-0.0802-0.0886-0.0980-0.1083-0.1197Columns73through80-0.1323-0.1462-0.1615-0.1785-0.1973-0.2180-0.2409-0.2663Columns81through88-0.2942-0.3251-0.3593-0.3971-0.4388-0.4849-0.5358-0.5921Columns89through96-0.6543-0.7230-0.7989-0.8828-0.9755-1.0780-1.1912-1.3163Columns97through101-1.4545-1.6073-1.7760-1.9626-2.1686改良欧拉法的精度比向前欧拉法更高。2-3分析应用题用以下三种不同的方法求下述微分方程的数值解，取画出解的图形，与准确值比拟并进展分析。1〕欧拉法；2〕改良欧拉法；2-4分析应用题考虑一个涉及到社会上与众不同的人的繁衍问题模型。假设在时刻(单位为年)，社会上有人口人，又假设所有与众不同的人与别的与众不同的人结婚后所生后代也是与众不同的人。而固定比例为的所有其他的后代也是与众不同的人。如果对所有人来说出生率假定为常数，又如果普通的人和与众不同的人的婚配是任意的，则此问题可以用微分方程表示为：其中变量表示在时刻社会上与众不同的人的比例，表示在时刻人口中与众不同的人的数量。1〕假定和，当步长为年时，求从到解的近似值，并作出近似解的曲线图形。2〕准确求出微分方程的解，并将你当时在分题(b)中得到的结果与此时的准确值进展比拟。1〕>>euler(0,50,50,0.01,inline('0.00