重庆大学
学生实验报告
实验课程名称数学实验
开课实验室DS1407
学院自动化年级2013 专业班自动化02班学生姓名侯刚学号********
开课时间2014 至2015 学年第二学期
数学与统计学院制
开课学院、实验室：数统学院DS1407实验时间：2014年4月3日
课程
名称数学实验实验项目
名称
种群数量的状态转移——
微分方程
实验项目类型
验证演示综合设计其他
指导
教师龚劬
成绩√
实验目的
[1] 归纳和学习求解常微分方程(组)的基本原理和方法；
[2] 掌握解析、数值解法，并学会用图形观察解的形态和进行解的定性分析；
[3] 熟悉MATLAB软件关于微分方程求解的各种命令；
[4] 通过范例学习建立微分方程方面的数学模型以及求解全过程；
基础实验
一、实验内容
1.微分方程及方程组的解析求解法；
2.微分方程及方程组的数值求解法——欧拉、欧拉改进算法；
3.直接使用MATLAB命令对微分方程(组)进行求解(包括解析解、数值解)；
4.利用图形对解的特征作定性分析；
5.建立微分方程方面的数学模型，并了解建立数学模型的全过程。
二、实验过程
1．求微分方程的解析解, 并画出它们的图形，
y’= y + 2x, y(0) = 1, 0<x<1;
（1）求解：
输入：dsolve('Dy=y+2*x','y(0)=1','x')
输出：ans=-2*x-2+3*exp(x)
（3）作图：
输入：>> x=0:0.1:1;
>> y2=-2*x-2+3*exp(x);
>> plot(x,y2)
输出：
图表 1 方程特解图形
分析：注意dsolve的用法。
2．用向前欧拉公式和改进的欧拉公式求方程y’= y - 2x/y, y(0) = 1 (0≤x≤1,h = 0.1) 的数值解，要求编写程序，并比较两种方法的计算结果，说明了什么问题？
（1）求解析解
输入： dsolve('Dy=y-2*x/y','y(0)=1','x')
输出： ans =(2*x+1)^(1/2)
（2）用向前欧拉公式和改进的欧拉公式求方程的数值解并与解析解作图比较
程序：
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));
k1=y2(k)-2*x1(k)/y2(k);
k2=y2(k)+h*k1-2*x1(k+1)/(y2(k)+h*k1);
y2(k+1)=y2(k)+h*(k1+k2)/2;
end
x1,y1,y2
x=0:0.1:1;
y=(2*x+1).^(1/2);
plot(x,y,x,y1,'o',x,y2,'+')
结果：
x1 =0 0.1000 0.2000 0.3000 0.4000 0.5000 0.6000 0.7000 0.8000 0.9000 1.0000
y1 =1.0000 1.1000 1.1918 1.2774 1.3582 1.4351 1.5090 1.5803
1.6498 1.7178 1.7848
y2 =1.0000 1.0959 1.1841 1.2662 1.3434 1.4164 1.4860 1.5525
1.6165 1.6782 1.7379
图表 2 向前欧拉公式和改进的欧拉公式所求方程数值解与解析解的比较
由图可得，改进后的欧拉公式求得的数值解更贴合解析解。
分析：注意向前欧拉与改进后的欧拉公式的不同。
3．Rossler 微分方程组：
当固定参数b=2, c=4时，试讨论随参数a 由小到大变化（如a ∈(0,0.65))而方程解的变化情况。
程序：
rossler.m:
function xdot=rossler(t,x)
xdot=[0,-1,-1;1,0.1,0;x(3),0,-4]*x+[0,0,2]'; fangchengzu.m: x0=[0 0 0.1];
[t,x]=ode45('rossler',[0,10],x0);
plot(t,x(:,1),'-',t,x(:,2),'.',t,x(:,3),'+') pause
plot3(x(:,1),x(:,2),x(:,3)) grid on 结果： a=0.1时：
a=0.25时：
⎪⎩
⎪
⎨⎧-+=+=--=)('''c x z b z ay
x y z y x
a=0.5时：
a=0.6时：
上述图形表示了a 由小到大变化时方程解的变化。
分析：注意xdot 的书写以及ode45的运用。
4.Apollo 卫星的运动轨迹的绘制
程序：
apollo.m:
function yp=apollo(t,x) u=1/82.45; u1=1-u;
r1=sqrt((x(1)+u)^2+x(3)^2);
11331
2
13
3
1212222
121()()
2,
2,
1/82.45,1,
(),()(0) 1.2,(0)0,(0)0,(0) 1.04935751
x x x y x r r y
y
y x y r r r x y r x y x x y y μμμμμμμμμμμ+-=+-
-
=-+-
-
==-=++=-+====-
r2=sqrt((x(1)-u1)^2+x(3)^2);
yp=[x(2);2*x(4)+x(1)-u1*(x(1)+u)/r1^3-u*(x(1)-u1)/r2^3;x(4);-2*x(2)+x(3)-u1*x(3)/r1^3-u*x(3)/r2^3];
weixing.m:
x0=[1.2;0;0;-1.04935751];
[t,x]=ode45('apollo',[0,20],x0);
plot(x(:,1),x(:,3))
xlabel('x')
ylabel('y')
title('Apollo卫星运动轨迹')
结果：
图表3 apollo卫星轨迹图
分析：注意求数值解时，高阶微分方程必须等价的变为一阶微分方程组。
应用实验（或综合实验）
一、实验内容
盐水的混合问题
一个圆柱形的容器，内装350升的均匀混合的盐水溶液。
如果纯水以每秒14升的速度从容器顶部流入，同时，容器内的混合的盐水以每秒10.5升的速度从容器底部流出。
开始时，容器内盐的含量为7千克。
求经过时间t后容器内盐的含量。
二、问题分析
（1）已知：水的密度为1kg/L,盐溶解度为36g。
可计算出7kg盐所需要的溶剂为194L水。
因此，由混合液体积即可知开始时刻的7kg盐是完全溶于水中的,并且没有饱和。
所以，整个过程为食盐水被再次稀释的过程，则不会出现有盐析出现象。
（2）由于容器的容积相对于单位时间内水的体积变化来说很大，所以可以忽略溶质盐在在不同浓度的水内扩散至均匀的时间。
根据在每个微小的时间段内，减少的盐加上容器内剩余的盐等于开始的盐量建立方程。
三、数学模型的建立与求解(一般应包括模型、求解步骤或思路，程序放在后面的附录中)
假设：1）温度对盐在水中的溶解度变化影响不大。
2）任意时刻容器内混合的、流出的盐水都均匀。
3）水流入及盐水流出的速度均为匀速。
设注水时间为t，t时刻时容器内含盐量为P(t)、容器内混合盐水的体积为V(t)，纯水流入容器的速度为v1,混合液流出的速度为v2。
可列出方程组：
P(t+△t)=P(t)-P(t)*v2*△t/V(t)
V(t)=V(t0)+(v1-v2)*t
V(t0)=350,P(0)=7,v1=14,v2=10.5
方程可化为：
dP/dt=-10.5*P(t)/(350+3.5*t),P(0)=7
用MATLAB求解该方程并作图。
四、实验结果及分析
求得方程的解析解为：
P(t)= 7000000/(t + 100)^3
曲线图像为：
图表 4 经过时间t后容器内盐的含量
五、附录（程序等）
y=dsolve('Dy=-14*y/(350+3.5*t)','y(0)=7','t')
ezplot('7000000/(t + 100)^3',[0,100])
xlabel('t')
ylabel('P(t)')
总结与体会
通过该实验的学习，掌握微分方程(组)求解方法（解析法、欧拉法、梯度法、改进欧拉法等），对常微分方程的数值解法有一个初步了解，同时学会使用MATLAB软件求解微分方程的基本命令，学会建立微分方程方面的数学模型。
教师签名
年月日。