实验五 种群数量的状态转移 - 微分方程（数学建模）

实验五 种群数量的状态转移——微分方程

一、实验目的及意义

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

[2] 掌握解析、数值解法，并学会用图形观察解的形态和进行解的定性分析； [3] 熟悉MATLAB软件关于微分方程求解的各种命令；

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

通过该实验的学习，使学生掌握微分方程(组)求解方法（解析法、欧拉法、梯度法、改进欧拉法等），对常微分方程的数值解法有一个初步了解，同时学会使用MATLAB软件求解微分方程的基本命令，学会建立微分方程方面的数学模型。这对于学生深入理解微分、积分的数学概念，掌握数学的分析思维方法，熟悉处理大量的工程计算问题的方法是十分必要的。

二、实验内容

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

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

3． 直接使用MATLAB命令对微分方程(组)进行求解(包括解析解、数值解)； 4． 利用图形对解的特征作定性分析；

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

三、实验步骤

1．开启软件平台——MATLAB，开启MATLAB编辑窗口； 2．根据微分方程求解步骤编写M文件 3．保存文件并运行；

4．观察运行结果(数值或图形)；

5．根据观察到的结果和体会写出实验报告。

四、实验要求与任务

根据实验内容和步骤，完成以下实验，要求写出实验报告（实验目的→问题→数学模型→算法与编程→计算结果→分析、检验和结论）

基础实验

1．求微分方程的解析解, 并画出它们的图形， y’= y + 2x, y(0) = 1, 0

y= dsolve('Dy-y-2*x=0','y(0)=1','x') 解出方程的解析解为 y =3*exp(x) - 2*x – 2

画出图形

X=linspace(0,1,100);

plot(x,3*exp(x) - 2*x - 2)

2．用向前欧拉公式和改进的欧拉公式求方程y’= y - 2x/y, y(0) = 1 (0≤x≤1,h = 0.1) 的数值解，要求编写程序，并比较两种方法的计算结果，说明了什么问题？ dsolve('Dy=y-2*x/y','y(0)=1','x') 求得解析解为 y=(2*x+1).^(1/2) 建立M文件weifen.m：

x1(1)=0;y1(1)=1;y2(1)=1;h=0.1; for n=1:10

x1(n+1)=x1(n)+h;

y1(n+1)=(1+h)*y1(n)-2*h*x1(n)/y1(n); k1=y2(n)-2*x1(n)/y2(n);

k2=y2(n)+h*k1-2*x1(n+1)/(y2(n)+h*k1); y2(n+1)=y2(n)+0.5*h*(k1+k2); end x1;y1; x=0:0.1:1; y=(2*x+1).^(1/2);

plot(x,y,x1,y1,'r',x1,y2,'k:')

由所作图形可明显看出用改进欧拉公式所得的微分曲线明显更接近微分方程的解析解曲线，所以用改进欧拉方程求得的数值解更精确

3．Rossler微分方程组：

?x'??y?z??y'?x?ay?z'?b?z(x?c)?

当固定参数b=2, c=4时，试讨论随参数a由小到大变化（如a∈(0,0.65))而方程解的变化情况，并且画出空间曲线图形，观察空间曲线是否形成混沌状？

function xdot=lorenz(t,x) a=0;0.01;0.65;

xdot=[-x(2)-x(3);x(1)+a*x(2);2+x(3)*(x(1)-4)]; end

x0=[0 0 0.1]';

[t,x]=ode45('lorenz',[0,10],x0);

plot(t,x(:,1),'-',t,x(:,2),'*',t,x(:,3),'+') pause

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

图中，x1的图形为实线(蓝），x2的图形为“*”线（绿）， x3的图形为“+”线（红）.取[t0，tf]=[0，10]。 若取[t0，tf]=[0，100]。图像为

4.Apollo卫星的运动轨迹的绘制

?(x??)?(x??1)

x?2y?x?13?,3rr 12

y??2x?y??1yr13??yr23,??1/82.45,?1?1??,r1?(x??)2?y2,r2?(x??1)2?y2x(0)?1.2,x(0)?0,y(0)?0,y(0)??1.04935751

令x(1)=x,x(2)=x’,x(3)=y,x(4)=y’ x’’=x(2)’,y’’=x(3)’ function f = weixing( t,x)

u=1/82.45;u1=1-u;r1=sqrt((x(1)+u)^2+x(3)^2);r2=sqrt((x(1)-u1)^2+x(3)^2);

f=[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)]; end

x0=[1.2,0,0,-1.04935751];[t,x]=ode45('weixing',[0,50],x0);plot(x(:,1),x(:,3))

应用实验

盐水的混合问题

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

记

开始的溶液体积为V0=350 L 在t时刻溶液体积为V(t) 纯水流入速度为

r= 14l/s 盐水流出速度为 k= 10.5 l/s

开始盐含量 Q0= 7kg 在t时刻含盐量为 Q(t) 在t时刻含浓度为 C(t)

在dt 时间内，流出容器的盐分dQ(t)为可认为在dt时间段内浓度不变为C(t) 有 dQ(t)=kC(t)dt 则 Q(t)=Q0-?kC(t)dt

两边求导得 Q(t)’=-kC(t)

由于 C(t)=Q(t)/V(t) V(t)=V0+(r 则 Q(t)’=-kQ(t)/[V0+(r?k)t

?k)t]

Matble 求解：

dsolve('Dq=-10.5*q/(350+3.5*t)','q(0)=7','t') ans =

7000000/(t + 100)^3

7?106可得 Q(t)?

(t?100)3t=linspace(0,500,10000);plot(t,7000000./(t + 100).^3)

画出Q(t)随时间变化曲线为

本文来源：https://www.bwwdw.com/article/k0ff.html