数学建模实验.docx

1、数学建模实验数学建模课程实验报告专题实验7班级数财系1班学号2011040123丛文实验题目常微分方程数值解实验目的1掌握用MATLAB求微分方程初值问题数值解的方法；2通过实例学习微分方程模型解决简化的实际问题；3了解欧拉方法和龙格库塔方法的基本思想。实验容(包括分析过程、方法、和代码，结果)1. 用欧拉方法和龙格库塔方法求下列微分方程初值问题的数值解，画出解的图形，对结果进行分析比较解；M文件function f=f(x,y)f=y+2*x;程序；clc;clear;a=0;b=1; %求解区间x1,y_r=ode45(f,a b,1); %调用龙格库塔求解函数求解数值解；% 以下利用Eu

2、ler方法求解y(1)=1;N=100;h=(b-a)/N;x=a:h:b;for i=1:N y(i+1)=y(i)+h*f(x(i),y(i);endfigure(1) plot(x1,y_r,r*,x,y,b+,x,3*exp(x)-2*x-2,k-);%数值解与真解图title(数值解与真解图);legend(RK4,Euler,真解);xlabel(x);ylabel(y);figure(2)plot(x1,abs(y_r-(3*exp(x1)-2*x1-2),k-);%龙格库塔方法的误差title(龙格库塔方法的误差)xlabel(x);ylabel(Error);figure(3

3、)plot(x,abs(y-(3*exp(x)-2*x-2),r-)%Euler方法的误差title(Euler方法的误差)xlabel(x);ylabel(Error);4. 单摆运动是一个我们熟悉的物理模型，可以看作工程技术中一些振动问题的简化，图8中一根长l的（无弹性的）细线，一端固定，另一端悬挂一质量为m的小球，在重力作用下小球处于竖直的平衡位置，使小球偏离平衡位置一根小的角度，然后让它无初速度的放开，小球就会沿圆弧摆动，在不考虑空气阻力的情况下建立关于时间t的微分方程，设l=25cm，在等于和两种情况下求方程的数值解，并与近似解比较。解；实验原理与数学模型：在小球摆动过程中的任一位置

4、，小球所受重力沿运动轨迹方向的分力为-mgsin（负号表示力的方向与的正方向相反），利用牛顿第二定律即得微分方程描述单摆运动规律的微分方程（1）是2阶微分方程，无解析解，但可用Matlab或其它软件编程求其数值解，但都需要先将它化成方程组的形式。令 则微分方程 (1) 化为 初始条件转化为 在前面的两式中，g=9.8，l=0.25，x10为10o=0.1745(弧度)及30o=0.5236(弧度)两种情况. 周期 根据上原理，可以建立模型，用matlab编程求解，过程如下，程序 以单摆的两个周期来计算作图。M文件；function dx=danbai(t,x)g=9.8;l=0.25;dx=x

5、(2);-g/l*sin(x(1);程序；一.当时ts=0:0.05:2;a0=0.1745;x0=a0,0;t,x=ode23(danbai,ts,x0);y=a0*cos(sqrt(40).*t);t,x(:,1),ysubplot(1,2,2),plot (t,x(:,1),-k*),title(摆角10度数值解)subplot(1,2,1), plot(t,y,b*),plot(t,y,-r*),title(摆角10度近似解)ans = 0 0.1745 0.1745 0.0500 0.1661 0.1658 0.1000 0.1415 0.1407 0.1500 0.1033 0.1

6、017 0.2000 0.0550 0.0525 0.2500 0.0015 -0.0018 0.3000 -0.0522 -0.0560 0.3500 -0.1009 -0.1046 0.4000 -0.1397 -0.1429 0.4500 -0.1649 -0.1669 0.5000 -0.1742 -0.1745 0.5500 -0.1667 -0.1647 0.6000 -0.1430 -0.1386 0.6500 -0.1054 -0.0987 0.7000 -0.0577 -0.0491 0.7500 -0.0044 0.0054 0.8000 0.0494 0.0594 0.8

7、500 0.0983 0.1075 0.9000 0.1377 0.1449 0.9500 0.1637 0.16791.0000 0.1738 0.1744二. 当时，ts=0:0.05:2;a0=0.5236;x0=a0,0;t,x=ode23(danbai,ts,x0);y=a0*cos(sqrt(40).*t);t,x(:,1),ysubplot(1,2,2),plot (t,x(:,1),-k*),title(摆角30度数值解)subplot(1,2,1), plot(t,y,b*),plot(t,y,-r*),title(摆角30度近似解)ans = 0 0.5236 0.5236

8、 0.0500 0.4993 0.4976 0.1000 0.4283 0.4223 0.1500 0.3170 0.3051 0.2000 0.1753 0.1577 0.2500 0.0169 -0.0054 0.3000 -0.1431 -0.1680 0.3500 -0.2893 -0.3139 0.4000 -0.4076 -0.4286 0.4500 -0.4874 -0.5009 0.5000 -0.5217 -0.5235 0.5500 -0.5075 -0.4942 0.6000 -0.4460 -0.4158 0.6500 -0.3425 -0.2963 0.7000 -0

9、.2064 -0.1473 0.7500 -0.0505 0.0162 0.8000 0.1102 0.1782 0.8500 0.2602 0.3225 0.9000 0.3852 0.4348 0.9500 0.4736 0.5039 1.0000 0.5177 0.5232 1.0500 0.5136 0.4905 1.1000 0.4617 0.4092 1.1500 0.3665 0.2873 1.2000 0.2365 0.1369 1.2500 0.0838 -0.0271 1.3000 -0.0770 -0.1883 1.3500 -0.2303 -0.3309 1.4000

10、-0.3613 -0.4407 1.4500 -0.4580 -0.5068 1.5000 -0.5116 -0.5226 1.5500 -0.5177 -0.4866 1.6000 -0.4755 -0.4023 1.6500 -0.3889 -0.2782 1.7000 -0.2654 -0.1264 1.7500 -0.1164 0.0379 1.8000 0.0437 0.1984 1.8500 0.1995 0.3393 1.9000 0.3361 0.4465 1.9500 0.4405 0.5094 2.0000 0.5036 0.5218从数据可以看出，角度为10o时精确（数值）解与近似解相差不大，而初始角度为30o时，随着时间的增加差别很大实验结果分析或者实验总结和体会龙格库塔方法和Euler方法求解常微分方程都能获得比较好的数值解，相比较而言龙格库塔方法的数值解的精度远远要比Euler方法的数值解的精度高。