计算方法实验七

浙江大学城市学院实验报告课程名称计算方法实验项目名称常微分方程初值问题的数值解法实验成绩指导老师（签名）日期2011-12-9一.实验目的和要求1．用Matlab软件掌握求微分方程数值解的欧拉方法；2．通过实例学习用微分方程模型解决简化的实际问题。二.实验内容和原理编程题2-1要求写出Matlab源程序(m文件)，并有适当的注释语句；分析应用题2-2，2-3，2-4，2-5要求将问题的分析过程、Matlab源程序和运行结果和结果的解释、算法的分析写在实验报告上。2-1编程编写用向前欧拉公式和改进欧拉公式求微分方程数值解的Matlab程序，问题如下：在区间,ab内(1)N个等距点处，逼近下列初值问题的解，并对程序的每一句添上注释语句。0(,)()yfxyaxbyay　　　　　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;fori=1:n;y(i+1)=y(i)+h*f(x(i),y(i));endplot(x,y)holdon%求微分方程的精确解x1=linspace(a,b,100);'精确解为's=dsolve(f1,b1,'x')symsxy1=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=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)holdon%求微分方程的精确解x1=linspace(a,b,100);'精确解为's=dsolve(f1,b1,'x')symsxy1=zeros(1,100);fori=1:100y1(i)=subs(s,x,x1(i));endplot(x1,y1,'r')title('红色代表精确解')2-2分析应用题假设等分区间数100n，用欧拉法和改进欧拉法在区间[0,10]t内求解初值问题()()20(0)10ytyty并作出解的曲线图形，同时将方程的解析解也画在同一张图上，并作比较，分析这两种方法的精度。(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*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（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*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分析应用题考虑一个涉及到社会上与众不同的人的繁衍问题模型。假设在时刻t(单位为年)，社会上有人口()xt人，又假设所有与众不同的人与别的与众不同的人结婚后所生后代也是与众不同的人。而固定比例为r的所有其他的后代也是与众不同的人。如果对所有人来说出生率假定为常数b，又如果普通的人和与众不同的人的婚配是任意的，则此问题可以用微分方程表示为：()(1())dptrbptdt其中变量()()()iptxtxt表示在时刻t社会上与众不同的人的比例，()ixt表示在时刻t人口中与众不同的人的数量。1）假定(0)0.01,0.02pb和0.1r，当步长为1h年时，求从0t到50t解()pt的近似值，并作出近似解的曲线图形。2）精确求出微分方程的解()pt，并将你当50t时在分题(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=Columns1through80.01000.01200.01400.01590.01790.01990.02180.0238Columns9through160.02570.02770.02960.03160.03350.03540.03740.0393Columns17through240.04120.04310.04500.04700.04890.05080.05270.0546Columns25through320.05640.05830.06020.06210.06400.06580.06770.0696Columns33through400.07140.07330.07510.07700.07880.08070.08250.0844Columns41through480.08620.08800.08980.09170.09350.09530.09710.0989Columns49through510.10070.10250.10432）dsolve('Dp=0.002-0.002*p','p(0)=0.01','t')ans=1-99/(100*exp(t/500))1-99/(100*exp(0.1))ans=0.1042与欧拉法求得的精确值0.1043差0.00012-4分析应用题设为小球作单摆运动时挂线与垂直方向的夹角，m为小球的质量，挂线长度为l，由牛顿第二定律即得微分方程sinmlmg其中表示求二阶导数，设小球初始偏离角度为0，且无初速，则方程的初始条件为0(0),(0)0，当0较大时，上述方程没有解析解。试用4阶Runge-Kutta公式（见下面matlab相关函数介绍）在0等于10和30两种情况下求解(设25lcm)，9.8g，并画出()t的图形。提示：先化为微分方程组，即令12,，可将微分方程化为微分方程组l1221sinddtdgdtl然后通过求解微分方程组得到1，即关于t的表达式。2-5分析应用题考虑固定端点横梁偏差的例子，假设服从于均匀负载。描述该物理现象的边值问题为22()02dSqxxlxldxEIEI其边界条件为(0)0和()0l。假设横梁特征如下：长度120lin，均匀负载强度q=100lb/ft，弹性系数72E=3.010lb/in.，端点受力S=1000lb，惯性的中心力矩4I=625in.。a.近似求每6in.长度时横梁的弯曲()x。b.实际的关系由12()()axaxxcecebxlxc给出，其中51c=9.24510450695458710，52c=9.50489549304541310，-4a=2.30940110，0.05b，1875000c。区间上的最大误差是否在0.2in之内？c.平衡定律要求0max()1/3000.0033xlx。横梁满足平衡吗？【MATLAB相关函数】求微分方程的解析解及其数值的代入dsolve(‘egn1’，‘egn2’，‘x’)subs(expr,{x,y,…},{x1,y1,…})其中‘egni’表示第i个方程，‘x’表示微分