计算方法实验七

1、浙江大学城市学院实验报告课程名称 计算方法 实验项目名称 常微分方程初值问题的数值解法 实验成绩 指导老师（签名 ） 日期 2011-12-9 一. 实验目的和要求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)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);endplot(x,y)hold 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); endplot(x1,y1,&#

3、39;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; 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)hold on% 求微分方程的精确解x1=linspace(a,b,100);'精确解为's=dsolve(f1,b1,'x')syms xy1

4、=zeros(1,100);for i=1:100 y1(i)=subs(s,x,x1(i);endplot(x1,y1,'r')title('红色代表精确解')2-2 分析应用题假设等分区间数，用欧拉法和改进欧拉法在区间内求解初值问题并作出解的曲线图形，同时将方程的解析解也画在同一张图上，并作比较，分析这两种方法的精度。(1)向前欧拉法>> euler(0,10,100,10,inline('y-20','x','y'),'Dy=y-20','y(0)=10')ans

5、=精确解为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 Columns 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 t

6、hrough 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 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

7、-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.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

8、 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('y-20','x','

9、;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 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

10、.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.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 Colu

11、mns 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.0294 -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

12、-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 -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 分析应用题

13、考虑一个涉及到社会上与众不同的人的繁衍问题模型。假设在时刻(单位为年)，社会上有人口人，又假设所有与众不同的人与别的与众不同的人结婚后所生后代也是与众不同的人。而固定比例为的所有其他的后代也是与众不同的人。如果对所有人来说出生率假定为常数，又如果普通的人和与众不同的人的婚配是任意的，则此问题可以用微分方程表示为： 其中变量表示在时刻社会上与众不同的人的比例，表示在时刻人口中与众不同的人的数量。 1）假定和，当步长为年时，求从到解的近似值，并作出近似解的曲线图形。2）精确求出微分方程的解，并将你当时在分题(b)中得到的结果与此时的精确值进行比较。1）>> euler(0,50,50,

14、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 through 8 0.0100 0.0120 0.0140 0.0159 0.0179 0.0199 0.0218 0.0238 Columns 9 through 16 0.0257 0.0277 0.0296 0.0316 0.0335 0.0354 0.0374 0.0

15、393 Columns 17 through 24 0.0412 0.0431 0.0450 0.0470 0.0489 0.0508 0.0527 0.0546 Columns 25 through 32 0.0564 0.0583 0.0602 0.0621 0.0640 0.0658 0.0677 0.0696 Columns 33 through 40 0.0714 0.0733 0.0751 0.0770 0.0788 0.0807 0.0825 0.0844 Columns 41 through 48 0.0862 0.0880 0.0898 0.0917 0.0935 0.095

16、3 0.0971 0.0989 Columns 49 through 510.1007 0.1025 0.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 分析应用题设为小球作单摆运动时挂线与垂直方向的夹角，为小球的质量，挂线长度为，由牛顿第二定律即得微分方程 其中表示求二阶导数，设小球初始偏离角度

17、为，且无初速，则方程的初始条件为，当较大时，上述方程没有解析解。试用4阶Runge-Kutta公式（见下面matlab相关函数介绍）在等于和两种情况下求解(设) ，并画出的图形。提示：先化为微分方程组，即令，可将微分方程化为微分方程组 然后通过求解微分方程组得到，即关于的表达式。2-5 分析应用题考虑固定端点横梁偏差的例子，假设服从于均匀负载。描述该物理现象的边值问题为 其边界条件为和。假设横梁特征如下：长度，均匀负载强度，弹性系数，端点受力，惯性的中心力矩。a. 近似求每6 in.长度时横梁的弯曲。b. 实际的关系由 给出，其中，。区间上的最大误差是否在0.2 in之内？c. 平衡定律要求。

18、横梁满足平衡吗？【MATLAB相关函数】n 求微分方程的解析解及其数值的代入 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 x>> subs('x2',x,2)ans = 4 s=dsolve(, , )ans = &

19、gt;> syms x>> subs(s,x,2)ans = -0.3721n 右端函数的自动生成 f= 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为函数的初值