重庆大学数学模型数学实验作业四讲解

开课学院、实验室：数统学院实验时间：2015年10月28日

课程名称数学实验实验项目

名称

种群数量的状态转移——

微分方程

实验项目类型

验证演示综合设计其他

指导

教师

肖剑成绩

实验目的

[1] 归纳和学习求解常微分方程(组)的基本原理和方法；

[2] 掌握解析、数值解法，并学会用图形观察解的形态和进行解的定性分析；

[3] 熟悉MATLAB软件关于微分方程求解的各种命令；

[4] 通过范例学习建立微分方程方面的数学模型以及求解全过程；

通过该实验的学习，使学生掌握微分方程(组)求解方法（解析法、欧拉法、梯度法、改进欧拉法等），对常微分方程的数值解法有一个初步了解，同时学会使用MATLAB软件求解微分方程的基本命令，学会建

立微分方程方面的数学模型。这对于学生深入理解微分、积分的数学概念，掌握数学的分析思维方法，熟

悉处理大量的工程计算问题的方法是十分必要的。

实验内容

1．微分方程及方程组的解析求解法；

2．微分方程及方程组的数值求解法——欧拉、欧拉改进算法；

3．直接使用MATLAB命令对微分方程(组)进行求解(包括解析解、数值解)；

4．利用图形对解的特征作定性分析；

5．建立微分方程方面的数学模型，并了解建立数学模型的全过程。

基础实验

一、问题重述

1．求微分方程的解析解, 并画出它们的图形，

y’= y + 2x, y(0) = 1, 0

2．用向前欧拉公式和改进的欧拉公式求方程y’= y - 2x/y, y(0) = 1 (0≤x≤1,h = 0.1) 的数值解，要求编写程序，并比较两种方法的计算结果，说明了什么问题？

3．Rossler微分方程组：

当固定参数

b=2, c=4

时，试讨论随参数a 由小到大变化（如

a ∈(0,0.65))

而方程解的变化情况，

并且画出空间曲线图形，观察空间曲线是否形成混沌状？

4.Apollo

卫星的运动轨迹的绘制

二、实验过程

1．编辑程序代码Untitle1:

s=dsolve('Dy=y+2*x'

,'y(0)=1'

,'x'

)

ezplot(s,[0,1])

运行结果如下：s =3*exp(x) - 2*x – 2

图形为：

1

1

331

2

13312

12

22

2

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 )

(''

'c x

z b

z ay x y z y x

2．编写程序代码Untitle2：

clc

y=dsolve('Dy=y-2*x/y','y(0)=1','x')

ezplot(y,[0,1])

hold on

x=[];x(1)=0;y1(1)=1;y2(1)=1;

for n=1:10

x(n+1)=x(n)+0.1;

y1(n+1)=1.1*y1(n)-0.2*x(n)/y1(n);

k1=y2(n)-2*x(n)/y2(n);

k2=y2(n)+0.1*k1-2*x(n)/(y2(n)+0.1*k1);

y2(n+1)=y2(n)+0.05*(k1+k2);

end

plot(x,y1,'k:',x,y2,'k-.')

运行得到y =(2*x + 1)^(1/2)，这是解析解。图像如下：

此题中向前欧拉公式更逼近解析解，其实，提高精度，即n的取值，两种方式都可以无限逼近解析解。

3.首先编辑函数m文件rossler.m：

function eq=rossler(t,x)

global a b c

b=2;c=4;

eq=[0 -1 -1;1 a 0;x(3) 0 -c]*x+[0;0;b];

然后在命令行窗口输入全局变量，并对a赋值，当a=0.1时：

global a b c

a=0.1;

x0=[0;0;0];

[t,x]=ode45('rossler',[0,1000],x0);

plot(t,x(:,1),'r',t,x(:,2),'k',t,x(:,3),'g') pause

plot3(x(:,1),x(:,2),x(:,3))

grid on

得到的图形为：

当a=0.2，积分区间改为[0,100]时：

当当a=0.6，积分区间改为[0,50]时：

由此一系列图可知此空间图线是混沌的。

4.编写函数m文件apollo.m:

function ep=apollo(t,y)

syms y1y2y3y4

u=1/82.45;u1=1-u;

r1=((y(1)+u)^2+y(3)^2)^(1/2);

r2=((y(1)-u1)^2+y(3)^2)^(1/2);

ep(1)=y(2);

ep(2)=2*y(4)+y(1)-u1*(y(1)+u)/r1^3-u*(y(1)-u1)/r2^3;

ep(3)=y(4);

ep(4)=-2*y(2)+y(3)-u1*y(3)/r1^3-u*y(3)/r2^3;

ep=[ep(1);ep(2);ep(3);ep(4)];

运行程序代码Untitle4：

[t,y]=ode45('apollo',[0,20],[1.2 0 0 -1.04935751]); y1=y(:,1);

y2=y(:,2);

y3=y(:,3);

y4=y(:,4);

plot(t,y1,t,y3)

grid on

pause

plot(y1,y3)

grid on

得到y1-t,y3-t的图像

y1-y3（Apollo卫星的运动轨迹）的图像：

应用实验（或综合实验）

一、问题重述

盐水的混合问题

一个圆柱形的容器，内装350升的均匀混合的盐水溶液。如果纯水以每秒14升的速度从容器顶部流入，同时，容器内的混合的盐水以每秒10.5升的速度从容器底部流出。开始时，容器内盐的含量为7千克。求经过时间t后容器内盐的含量。

二、问题分析

由题意可以知道，此题中容器内的盐含量以及浓度随着时间在不停变化的，在流入到流出的过程

中，由于混合在水中的盐含量是不同的，所以溶解于水中的盐的量每一时刻都是不同的，流出的量随时间

也是不断变化的。可以选取一个无限小的时间微元进行讨论。

三、数学模型的建立与求解

假设在t时刻到t+△t（△t足够小）时刻时，由于时间变化非常微小，可以认为这个△t时间内，容器内溶液浓度没有发生变化浓度c(t)=c(t+△t)=m(t)/V(t),则这个过程中盐减少的质量为△m=m(t+△t)-m(t)=-v2*c(t)*△t，V(t)=V0+v1*t-v2*t。其中c(t)表示t时刻容器内盐的浓度；m(t)表示t时刻容器内盐的质量；V(t)表示t时刻容器内水的体积；v1,v2分别表示流入流出水的速度。

由于△t足够小，得微分方程m’(t)=-v2*m(t)/(V(0)+v1*t-v2*t)。带入数据得

dm/dt=3m/(100+t)

在MA TLAB中建立脚本Untitle.m

m=dsolve ('Dm=-3*m/(100+t)','m(0)=7','t')

ezplot(m,[0,150])

grid on

四、实验结果及分析

运行程序得到解析解m =7000000/(t + 100)^3，以及它的图像

可以看出t越大时y是接近于0的，但是不等于0，与真实情况相符。

教师签名

年月日