重庆大学数学实验3微分方程

1、重 庆 大 学学 生 实 验 报 告实验课程名称 数学实验 开课实验室 DS1422 学 院 学 生 姓 名 开 课 时 间 总 成 绩教师签名数学与统计 学 院 制开课学院、实验室： 数统学院 DS1422 实验时间 ：课程名称数学实验实验项目名 称实验项目类型验证演示综合设计其他指导教师肖剑成 绩实验目的1 归纳和学习求解常微分方程(组)的基本原理和方法；2 掌握解析、数值解法，并学会用图形观察解的形态和进行解的定性分析；3 熟悉MATLAB软件关于微分方程求解的各种命令；4 通过范例学习建立微分方程方面的数学模型以及求解全过程； 通过该实验的学习，使学生掌握微分方程(组)求解方法（解析法

2、、欧拉法、梯度法、改进欧拉法等），对常微分方程的数值解法有一个初步了解，同时学会使用MATLAB软件求解微分方程的基本命令，学会建立微分方程方面的数学模型。这对于学生深入理解微分、积分的数学概念，掌握数学的分析思维方法，熟悉处理大量的工程计算问题的方法是十分必要的。基础实验一、实验内容1 微分方程及方程组的解析求解法；2 微分方程及方程组的数值求解法欧拉、欧拉改进算法；3 直接使用MATLAB命令对微分方程(组)进行求解(包括解析解、数值解)；4 利用图形对解的特征作定性分析；5 建立微分方程方面的数学模型，并了解建立数学模型的全过程。二、实验过程（一般应包括实验原理或问题分析，算法设计、程序

3、、计算、图表等， 实验结果及分析） 1求微分方程的解析解, 并画出它们的图形， y= y + 2x, y(0) = 1, 0<x<1; y+ycos(x) = 0, y(0)=1, y(0)=0; （1）通解y=dsolve('Dy=y+2*x' ,'x') y = -2*x-2+exp(x)*C1特解y1=dsolve('Dy=y+2*x','y(0)=1','x') y1 = -2*x-2+3*exp(x)作图（y =3*exp(x)-2*x2）x=0:0.02:1;y1=3*exp(x)-2*x-

4、2;plot(x,y1,'b-')（2）算法设计通解y=dsolve('D2y+y*cos(x)=0') y = C19*exp(t*(-cos(x)(1/2) + C20/exp(t*(-cos(x)(1/2)特解 y=dsolve('D2y+y*cos(x)=0','y(0)=1,Dy(0)=0','x')结果 y = empty sym 分析：此微分方程没有解析解.该方程数值解的求解过程：设y1=y,y2=y 则可得y1=y2，y2+y1cos(x)=0首先建立M-文件：function f=fun(x,y)

5、f=y(2);-y(1)*cos(x);输入命令：x,y=ode23('fun',0,1,1,0)结果：x = 0 0.0001 0.0005 0.0025 0.0125 0.0625 0.1541 0.2541 0.3541 0.4541 0.5541 0.6541 0.7541 0.8541 0.9541 1.0000y = 1.0000 0 1.0000 -0.0001 1.0000 -0.0005 1.0000 -0.0025 0.9999 -0.0125 0.9980 -0.0624 0.9882 -0.1529 0.9680 -0.2487 0.9386 -0.33

6、97 0.9003 -0.4243 0.8540 -0.5011 0.8004 -0.5693 0.7404 -0.6280 0.6751 -0.6772 0.6053 -0.7168 0.5721 -0.7319作图：y1=y(:,1);>> y2=y(:,2);>> plot(x,y1,'b*',x,y2,'r'),gtext('y1'),gtext('y2')2用向前欧拉公式和改进的欧拉公式求方程y= y - 2x/y, y(0) = 1 (0x1,h = 0.1) 的数值解，要求编写程序，并比较两种

7、方法的计算结果，说明了什么问题？程序：向前欧拉和向后欧拉法:x1(1)=0;y1(1)=1;y2(1)=1;h=0.1; for k=1;10 x1(k+1)=x1(k)+h;y1(k+1)=y1(k)+h*(y1(k)-2*x1(k)/y1(k) y2(k+1)=y2(k)+h*(y2(k+1)-2*x1(k+1)/y2(k+1);endx1,y1,y2ans = 10y1 = 1.0000 1.1000 1.0000 1.0000 0.9999 0.9980 0.9882 0.9680 0.9386 0.9003 0.8540 0.8004 0.7404 0.6751 0.6053 0.5

8、721x1 = 0 0.1000y1 = 1.0000 1.1000 1.0000 1.0000 0.9999 0.9980 0.9882 0.9680 0.9386 0.9003 0.8540 0.8004 0.7404 0.6751 0.6053 0.5721y2 = 1.0000 251.0000 -0.0005 -0.0025 -0.0125 -0.0624 -0.1529 -0.2487 -0.3397 -0.4243 -0.5011 -0.5693 -0.6280 -0.6772 -0.7168 -0.7319改进欧拉公式法：算法设计：x1(1)=0;y1(1)=1;h=0.1;f

9、or k=1:10 x1(k+1)=x1(k)+h; y1(k+1)=y1(k)+1/2*h*(y1(k)-2*x1(k)/y1(k)+y1(k+1)-2*x1(k+1)/y1(k+1);endx1,y1结果：x1 = Columns 1 through 6 0 0.1000 0.2000 0.3000 0.4000 0.5000 Columns 7 through 11 0.6000 0.7000 0.8000 0.9000 1.0000y1 = 1.0000 1.0957 1.1828 1.2624 1.3345 1.3989 1.4544 1.4994 1.5313 1.5462 1.5

10、380 0.8004 0.7404 0.6751 0.6053 0.57213Rossler微分方程组;当固定参数b=2, c=4时，试讨论随参数a由小到大变化（如a(0,0.65)而方程解的变化情况，并且画出空间曲线图形，观察空间曲线是否形成混沌状？建立M-file文件function f=fun(t,x)f=0,-1,-1;1,0.1,0;x(3),0,-4*x+0,0,2'输入命令：x0=0 0 0.1;t,x=ode45('fun',0,10,x0);plot(t,x(:,1),'r',t,x(:,2),'b-',t,x(:,3)

11、,'-')pauseplot3(x(:,1),x(:,2),x(:,3),grid on当a=0.2时 当a=0.3时 当a=0.4时当a= 0.5当a=0.6时结论：空间曲线不会形成混沌状4.Apollo卫星的运动轨迹的绘制则编写M文件函数：function f=fun(x,y)u=1/82.45,u1=1-u,r1=sqrt(y(1)+u)2+y(3)2),r2=sqrt(y(1)-u1)2+y(3)2)f=y(2);2*y(4)+y(1)-u1*(y(1)+u)/(r1)3-u*(y(1)-u1)/(r2)3;y(4);-2*y(2)+y(3)-u1*y(3)/(r1)3

12、-u*y(3)/(r2)3输入：x,y=ode45('fun',0,10,1.2,0,0,-1.04935751)plot(y(:,1),y(:,3)'),grid on作图：应用实验（或综合实验）一、实验内容5盐水的混合问题一个圆柱形的容器，内装350升的均匀混合的盐水溶液。如果纯水以每秒14升的速度从容器顶部流入，同时，容器内的混合的盐水以每秒10.5升的速度从容器底部流出。开始时，容器内盐的含量为7千克。求经过时间t后容器内盐的含量。二、问题分析 已知：水的密度为1kg/L,盐溶解度为36g(在假设（1）中)由题意可以知道，在流入到流出的过程中，由于混合在水中的盐

13、含量是不同的，所以溶解于水中的盐的量每一时刻都是不同的。而所求的这个瞬时量即为微分方程的解。可计算出7kg盐所需要的溶剂为194L水。因此，由混合液体积即可知开始时刻的7kg盐是完全溶于水中的,并且没有饱和。所以，整个过程为食盐水被再次稀释的过程，则不会出现有盐析出现象。三、数学模型的建立与求解(一般应包括模型、求解步骤或思路，程序放在后面的附录中) 假设：1）温度对盐在水中的溶解度变化影响不大， 2）任意时刻容器内混合的、流出的盐水都均匀3）水流入及盐水流出的速度均为匀速的。4）假设注水时间为t、t时刻时容器内含盐量为P(t)、纯水流入容器的速度为v1,混合液流出的速度为v2。计算过程如下：1）t时刻容器内含盐量可表示为P(t)2）纯水以v1=14L/s的速度流入容器，混合液以v2=10.5L/s的速度流出容器，则根据溶解度公式：，可推出7kg食盐是完全溶解在水中的3）根据假设1）、2），则可列出关于盐在水中的浓度的比例式，即为（假设其中足够小）4）根据上述比例式即可得到相应的方程为：而使可得微分方程：。5）求解该微分方程当时的曲线图像四、实验结果及分析 编写M文件，并根据得出的结果作图编写M文件：y=dsolve('Dy=-10.5*y/(350+3.5*t)','y(0)=7','t&#