数学建模作业微分方程实验北京工业大学.docx
《数学建模作业微分方程实验北京工业大学.docx》由会员分享，可在线阅读，更多相关《数学建模作业微分方程实验北京工业大学.docx（28页珍藏版）》请在冰点文库上搜索。
数学建模作业微分方程实验北京工业大学
2微分方程实验
1、微分方程稳定性分析
绘出下列自治系统相应的轨线，并标出随t增加的运动方向，确定平衡点，并按稳定的、渐近稳定的、或不稳定的进行分类：
解：
（1）由f（x）=x=0，f（y）=y=0；可得平衡点为（0,0），
系数矩阵
，求得特征值λ1=1，λ2=1；
p=-（λ1+λ2）=-2<0，q=λ1λ2=1>0；对照稳定性的情况表，可知平衡点（0，0）是不稳定的。
图形如下：
（2）如上题可求得平衡点为（0,0），特征值λ1=-1，λ2=2；
p=-（λ1+λ2）=-1<0，q=λ1λ2=-2<0；对照稳定性的情况表，可知平衡点（0，0）是不稳定的。
其图形如下：
（3）如上题可求得平衡点为（0,0），特征值λ1=0+1.4142i，λ2=0-1.4142i；
p=-（λ1+λ2）=0，q=λ1λ2=1.4142>0；对照稳定性的情况表，可知平衡点（0，0）是不稳定的。
其图形如下：
（4）如上题可求得平衡点为（1,0），特征值λ1=-1，λ2=-2；
p=-（λ1+λ2）=3>0，q=λ1λ2=2>0；对照稳定性的情况表，可知平衡点（1，0）是稳定的。
其图形如下：
2、种群增长模型
一个片子上的一群病菌趋向于繁殖成一个圆菌落.设病菌的数目为N，单位成员的增长率为r1，则由Malthus生长律有
，但是，处于周界表面的那些病菌由于寒冷而受到损伤，它们死亡的数量与N1/2成比例，其比例系数为r2，求N满足的微分方程.不用求解，图示其解族.方程是否有平衡解，如果有，是否为稳定的?
解：
由题意很容易列出N满足的微分方程：
令
=0，可求得方程的两个平衡点N1=0,N2=
进而求得
令
可求得N=
则N=N1，N=N2，N=
可以把第一象限划为三部分，且从下到上三部分中分别有
；
；
。
则可以画出N（t）的图形，即微分方程的解族，如下图所示：
由图形也可以看出，对于方程的两个平衡点，其中N1=0是不稳定的；N2=
是稳定的。
3、有限资源竞争模型
1926年Volterra提出了两个物种为共同的、有限的食物来源而竞争的模型
假设
，称
为物种i对食物不足的敏感度，
（1）证明当x1（t0）>0时，物种2最终要灭亡;
（2）用图形分析方法来说明物种2最终要灭亡.
解：
（1）由上述方程组f（x1）=
=0，
f（x2）=
=0，可得方程的平衡点为P0（0,0），P1（
0），P2（0,
）.
对平衡点P0（0,0），
系数矩阵
则p=-（b1+b2）<0，所以该平衡点不稳定。
对平衡点P1（
0），
系数矩阵
则p=
，q=
，
由题意
，x1（t0）>0，可以得出p>0,q>0,因此该平衡点是稳定的。
即
时，
，说明物种2最终要灭亡。
对平衡点P2（0,
），同理可以得到q<0，在该平衡点不稳定。
因此，在
，x1（t0）>0的条件下，物种2最终要灭亡。
（2）对于线性方程组
在平面上匹配两条直线l1和l2，由题意
，x1（t0）>0，可将第一象限分为三个区域。
在最左边区域，
都大于0；在中间区域，
都小于0，在最右边区域，
分别是大于0和小于0.，由各区域中
的取值可得到如下图形：
由图也可以看出，随着时间的增加，物种1最终能达到稳定值
，物种2最终要灭亡。
4、蝴蝶效应与混沌解
考虑Lorenz模型
其中σ=10，ρ=28，β=8/3，且初值为,x1（0）=x2（0）=0，x3（0）=ε，ε为一个小常数，假设ε=10-10，且0≤t≤100。
（1）用函数ode45求解，并画出x2~x1,x2~x3,x3~x1的平面图；
（2）适当地调整参数σ，ρ，β值，和初始值x1（0），x2（0）=0，x3（0），重复一的工作，看有什么现象发生。
解：
（1）编写Lorenz函数，
functionxdot=lorenz1（t,x,b,a,c）
xdot=[-b*x
（1）+x
（2）*x（3）;
-a*x
（2）+a*x（3）;
-x
（1）*x
（2）+c*x
（2）-x（3）];
对各参数赋值并用ode45函数求解，可得数值解：
Columns1through9
00.12500.25000.37500.50000.53520.57050.60570.6409
00.00000.0000-0.0000-0.00000.00000.00000.00000.0000
00.00000.0000-0.0000-0.00000.00000.00000.0000-0.0000
0.0000-0.0000-0.00000.00000.00000.0000-0.00000.00000.0000
Columns10through18
0.67610.71140.74660.78180.83080.87970.92860.97761.0105
0.00000.00000.00000.00000.00000.00000.00000.00000.0000
0.00000.00000.00000.00000.00000.00000.00000.00000.0000
0.00000.00000.00000.00000.00000.00000.00000.00000.0000
Columns19through27
1.04341.07631.10921.14211.17501.20791.24091.27971.3186
0.00000.00000.00000.00000.00000.00000.00000.00000.0000
0.00000.00000.00000.00000.00000.00000.00000.00010.0001
0.00000.00000.00000.00000.00000.00010.00010.00020.0002
Columns28through36
1.35751.39641.42461.45281.48101.50921.53741.56561.5938
0.00000.00000.00000.00000.00000.00000.00000.00000.0000
0.00020.00030.00040.00050.00080.00110.00150.00210.0029
0.00040.00060.00080.00120.00160.00230.00320.00450.0063
...............
Columns5590through5598
99.575199.592199.609099.626099.646299.666499.686799.706999.7338
16.945716.526116.201015.985415.897816.034816.447617.193318.7894
-3.3551-3.7119-4.1098-4.5568-5.1636-5.8601-6.6527-7.5438-8.8677
-5.3476-5.9274-6.5941-7.3519-8.3766-9.5370-10.8209-12.1944-14.0725
Columns5599through5607
99.760899.787899.814899.830699.846599.862399.878299.894099.9099
21.218624.470928.262930.535432.653334.464535.846636.714937.0528
-10.3044-11.7459-12.9932-13.5361-13.8800-13.9854-13.8351-13.4302-12.7919
-15.7415-16.8309-16.9274-16.3648-15.3302-13.8707-12.0805-10.1013-8.1005
Columns5608through5613
99.925799.941699.956299.970899.9854100.0000
36.895736.323935.524834.543533.448132.2971
-11.9606-10.9899-10.0237-9.0332-8.0563-7.1223
-6.2187-4.5596-3.2818-2.2590-1.4835-0.9275
x2~x1,x2~x3,x3~x1的平面图分别如下：
（2）令参数σ，ρ，β值各减1，初始值x1（0），x2（0）不变，x3（0）=10-8
分别得到得到x2~x1,x2~x3,x3~x1的平面图如下：
可以看出，参数σ，ρ，β值各减1，初始值x1（0），x2（0）不变，x3（0）数值变为x3（0）=10-8，参数和初始值很小的改变，就会导致最后图形发生很大的变化。
5、用微分方程考察共振现象
设物体沿x轴运动（如图所示）其平衡位置取为原点0，物体的质量为1，在时间t物体的位置为x（t）其所受的恢复为（如弹性力等）与物体所在位置的坐标成正比，即k2x,其中常数k称为恢复系数，运动过程所受的阻力（由于介质及摩擦等）设与速度成正比，即
，h>0，称为阻尼系数。
（1）根据Newton第二定律，建立相应的微分方程.不妨设初始位置为1，初始速度为0,取k=2,h=0（当h=0称为简谐振动的方程）和h=0.1，用Matlab软件得到相应的数值解，并在t-x平面上画出x（t）的图形。
（2）如果物体还受到附加外力的干扰，且外力是一个依据时间t的函数f（t）（设f（t）=Bsinwt），建立相应的微分方程（该方程称为强迫振动方程）.在上述参数不变的情况下，取振幅B=1，分别取w=1,1.2,1.4,1.6,1.8,2.0,2.2,2.4,2.6,2.8,3.0，用Matlab软件得到相应的数值解，并在t-x平面上画出x（t）的图形。
（3）分别对上述图形讲行分析，并解释为什么会出现这些现象。
解：
（1）根据Newton第二定律，F=ma，可得微分方程：
由题意知边界值：
x（0）=1,x’（0）=0,令y1=x，y2=
，
可将二阶方程转化为：
其中，m=1；g=9.8；k=2
当h=0时，由matlab解得数值解：
Columns1through9
00.00000.00000.00000.00000.00010.00010.00020.0002
1.00001.00001.00001.00001.00001.00001.00001.00001.0000
Columns10through18
0.00040.00060.00090.00110.00220.00320.00430.00540.0108
1.00001.00001.00001.00001.00001.00001.00011.00011.0003
Columns19through27
0.01620.02160.02710.05410.08120.10830.13530.21550.2956
1.00081.00141.00211.00851.01911.03391.05281.13261.2461
............
Columns100through108
8.23218.35938.46158.56368.66588.76798.87668.98529.0939
3.50193.21612.95052.66412.36882.07691.78341.52141.3034
Columns109through117
9.20259.31429.42609.53779.64949.73719.82479.912410.0000
1.13931.03451.00011.03831.14671.27731.44381.64121.8634
h=0时，x（t）图形如下所示：
当h=0.1时，由matlab解得数值解：
Columns1through9
00.00000.00000.00000.00000.00010.00010.00020.0002
1.00001.00001.00001.00001.00001.00001.00001.00001.0000
Columns10through18
0.00040.00060.00090.00110.00220.00320.00430.00540.0108
1.00001.00001.00001.00001.00001.00001.00011.00011.0003
Columns19through27
0.01620.02160.02710.05410.08120.10830.13530.21560.2958
1.00081.00141.00211.00851.01901.03361.05231.13081.2417
.............
Columns109through117
8.56318.66908.77498.88098.98689.12239.25799.39349.5289
2.58572.45632.32942.21062.10491.99521.92131.88801.8962
Columns118through121
9.64679.76459.882210.0000
1.93532.00162.09092.1979
h=0.1时，x（t）图形如下所示：
（2）如果还受到外力干扰，则微分方程变为：
将其化为一阶方程组：
其中，m=1；g=9.8；k=2；B=1；h=0.1；
用Matlab解得数值解：
w=1时
Columns1through9
00.00000.00000.00000.00000.00010.00010.00020.0002
1.00001.00001.00001.00001.00001.00001.00001.00001.0000
............
Columns118through121
9.82079.88049.940210.0000
1.89731.92101.95031.9846
w=1.2时
Columns1through9
00.00000.00000.00000.00000.00010.00010.00020.0002
1.00001.00001.00001.00001.00001.00001.00001.00001.0000
............
Columns118through121
9.82819.88549.942710.0000
1.69991.75681.81991.8886
w=1.4时
Columns1through9
00.00000.00000.00000.00000.00010.00010.00020.0002
1.00001.00001.00001.00001.00001.00001.00001.00001.0000
............
Columns118through125
9.62369.73099.83829.94549.95919.97279.986410.0000
2.24012.31712.40842.51032.52382.53732.55102.5647
w=1.4时
Columns1through9
00.00000.00000.00000.00000.00010.00010.00020.0002
1.00001.00001.00001.00001.00001.00001.00001.00001.0000
............
Columns118through121
9.78239.85499.927410.0000
2.11502.06762.02851.9979
w=1.6时
Columns1through9
00.00000.00000.00000.00000.00010.00010.00020.0002
1.00001.00001.00001.00001.00001.00001.00001.00001.0000
............
Columns118through121
9.69339.79559.897810.0000
0.80620.75030.75620.8231
w=1.8时
Columns1through9
00.00000.00000.00000.00000.00010.00010.00020.0002
1.00001.00001.00001.00001.00001.00001.00001.00001.0000
..............
Columns118through121
9.77529.85019.925110.0000
0.86291.08501.33831.6174
w=2.0时
Columns1through9
00.00000.00000.00000.00000.00010.00010.00020.0002
1.00001.00001.00001.00001.00001.00001.00001.00001.0000
.............
Columns118through121
9.75389.83599.917910.0000
2.32782.60232.87053.1243
.............
w=3.0时
Columns1through9
00.00000.00000.00000.00000.00010.00010.00020.0002
1.00001.00001.00001.00001.00001.00001.00001.00001.0000
...........
Columns118through121
9.74959.83309.916510.0000
2.24142.33032.41452.4915
如上图所示，w分别为1.0；1.2；...；3.0时x（t）的图形。
当h=0时，即无阻尼振动的情况下，得到t（x）图形如下所示：
（3）由上述图形可以看出，不考虑外加作用力时，当没有摩擦等阻力的时候，简谐运动一直进行下去，当有摩擦等阻力的时候，其振幅会逐渐的减小，最后趋近于零，但是周期不会发生变化。
在外加一个作用力之后，整个运动情况就发生了很大的变化。
首先刚开始的振幅和周期都变化很大，后来随着时间的推移，他们也逐渐开始趋于稳定。
但是，稳定之后的振幅和周期都各不相同，与外力的频率有直接的关系，这就是受迫振动，不仅跟系统本身的性质有关，还和外界干扰的情况有很大关系。
恢复系数k与w的值越接近时，随着时间增加，物体的振动振幅会一直增大，当w=k=2时，物体的振幅增大的速度最快，可以预见，在t趋于无穷大时，振幅也趋于无穷大，这就是共振现象。
也就是当外力与系统处于共振时，会引起振幅无限增大的振动，这在机械和建筑中一般是必须严格避免的。
6、人口模型与曲线拟合问题
表2.1列出的是美国1790一1980年的人口统计表.
（1）试用Malthus人口模型，按三段时间（1'790-1850,1850-1910,1910-1970）分别确定其增长率r。
并将数据和不同r值的三段曲线画在同一张图上.
（2）利用Logistic模型，重新确定固有增长率r和最大容量Nm，作图，再利用该模型得到的结果预测1990年的美国人口数。
解：
（1）分段研究，我们先求出增长率，编写命令如下：
t=10:
10:
70;
x=[3.95.37.29.612.917.123.2];
p=polyfit（t,log（x）,1）;p
得到结果p=0.0296
即1790——1850年时间段的增长率是0.0296
同理可以得到另外两段的增长率，1850-1910年为0.0226
1910——1970年为0.0129
每次画出图形，使用以下命令
plot（t,x,'bx'）;
holdon
在以上每一次求时顺便画出图形，得到最后的总图如下
（2）利用Logistics模型，将以上所有数据进行拟合，编写程序如下
t=0:
10:
190;
x=[3.95.37.29.612.917.123.231.438.650.262.976.092.0106.5123.2131.7150.7179.3204.0226.5];
f=inline（'p
（1）./（1+p
（2）*exp（-p（3）.*t））','p','t'）;
p=lsqcurvefit（f,[300,50,0.02],t,x）;
得到结果
人口最大容量Nm是360.3560（百万），增长率r是0.0234
如下图所示：
输入
t_pre=200;
x_pre=p
（1）./（1+p
（2）*exp（-p（3）*t_pre））;
可得最后的结果，即1990年的人口总数
x_pre=241.7704万人。
如下图：
7、加分实验（餐厅废物的堆肥优化问题）
一家环保餐厅用微生物将剩余的食物变成肥料。
餐厅每天将剩余的食物制成桨状物并与蔬菜下脚及少量纸片混合成原料，加入真菌菌种后放入容器内。
真菌消化这此混合原料，变成肥料，由于原料充足，肥料需求旺盛，餐厅希望增加肥料产量。
由于无力购置新设备，餐厅希望用增加真菌活力的办法来加速肥料生产.试通过分析以前肥料生产的记录（如表2.2所示），建立反映肥料生成机理的数学模型
升级会员 