当前位置: 首页 > >

数学建模作业、微分方程实验、北京工业大学

发布时间:

2 微分方程实验

1、微分方程稳定性分析 绘出下列自治系统相应的轨线,并标出随 t 增加的运动方向,确定*衡点, 并按稳定的、渐*稳定的、或不稳定的进行分类:
? dx ? dx ? dx ? dx ? x, ? ? x, ? y, ? ? x +1, ? dt ? dt ? dt ? dt ? ? ? ? (1) ? (2) ? (3) ? (4) ? dy dy dy ? ? ? ? dy ? ? 2 y. ? y; ? 2 y; ? ? 2 x; ? dt ? dt ? dt ? dt ? ? ? ?

解: (1)由 f(x)=x=0,f(y)=y=0;可得*衡点为(0,0) , 系数矩阵 A
?1 ? ? ?0 0? ? 1?

,求得特征值λ 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 生长律有
dN dt ? r1 ? N

,但是,处于周界表面的

那些病菌由于寒冷而受到损伤, 它们死亡的数量与 N1/2 成比例, 其比例系数为 r2, 求 N 满足的微分方程.不用求解,图示其解族.方程是否有*衡解,如果有,是否 为稳定的? 解:由题意很容易列出 N 满足的微分方程: 令
dN dt
1

? r1 N ? r2 N

2

? f (N )

f (N )

=0,可求得方程的两个*衡点 N1=0,N2= r2 2
d N dt
2 2

/ r1

2

进而求得
2

? ( r1 ?

1 2

?

1 2

1

r2 N

) ? ( r1 N ? r2 N

2

)



d N dt
2

? 0

可求得 N= r2 2

/ 4 r1

2

则 N=N1,N=N2,N= r2 2 别有
dN dt ? 0,

/ 4 r1

2

可以把第一象限划为三部分,且从下到上三部分中分
? 0,

d N dt
2

2

? 0, ;

dN dt

d N dt
2

2

? 0;

dN dt

? 0,

d N dt
2

2

? 0。

则可以画出 N(t)的图形,即微分方程的解族,如下图所示:

由图形也可以看出,对于方程的两个*衡点,其中 N1=0 是不稳定的;N2=
r2 / r1
2 2

是稳定的。

3、有限资源竞争模型 1926 年 Volterra 提出了两个物种为共同的、有限的食物来源而竞争的模型
? d x1 ? [ b1 ? ? 1 ( h1 x1 ? h 2 x 2 )] x1 ? dt ? ? ? d x 2 ? [ b ? ? ( h x ? h x )] x 2 2 1 1 2 2 2 ? dt ?

假设

b1

?1

?

b2

?2

,称

bi

?i

为物种 i 对食物不足的敏感度,

(1)证明当 x1(t0)>0 时,物种 2 最终要灭亡; (2)用图形分析方法来说明物种 2 最终要灭亡. 解: (1)由上述方程组 f(x1)= [ b1 ? ?1 ( h1 x1 ? h2 x 2 )] x1 =0, f(x2)= [ b2
b2

? ? 2 ( h1 x1 ? h 2 x 2 )] x 2 =0,可得方程的*衡点为

P0(0,0) 1( ,P

b1

? 1 h1

,0) ,

P2(0,

? 2 h2

).

对*衡点 P0(0,0) , 系数矩阵 A ? ?
? ? b1 ? h1 ? 1 x 1 ? h 2 ? 1 x 2 ? h1 ? 2 x 2 ? ? b1 ? ? ? b 2 ? h1 ? 2 x1 ? h1 ? 2 x 2 ? ? 0 ? h 2 ? 2 x1 0 ? ? b2 ?

则 p=-(b1+b2)<0,所以该*衡点不稳定。
b1

对*衡点 P1(

? 1 h1

,0) ,
? ?b ? 1 ? h 2 ? 2 x1 ? ? ? ? b 2 ? h1 ? 2 x1 ? h1 ? 2 x 2 ? ? ? 0 ?
)

系数矩阵 A

? b1 ? h1 ? 1 x1 ? h 2 ? 1 x 2 ? ? ? h1 ? 2 x 2 ?

h 2 b1 ? h1 ? ? b1 ? 2 ? b2 ? ? ?1 ? ?

则 p= b1 ? b 2
b1

?

b1 ? 2

?1

,q= ? b1 ( b 2

?

b1 ? 2

?1



由题意

?1

?

b2

?2

,x1(t0)>0,可以得出 p>0,q>0,因此该*衡点是稳定的。

即t

? ?

时, ( x1 ( t ), x 2 ( t )) ?

(

b1

? 1 h1

, 0)

,说明物种 2 最终要灭亡。

对*衡点 P2(0,

b2

? 2 h2

) ,同理可以得到 q<0,在该*衡点不稳定。

因此,在

b1

?1

?

b2

?2

,x1(t0)>0 的条件下,物种 2 最终要灭亡。
? b1 ? ? 1 ( h1 x 1 ? h 2 x 2 ) ? 0 ? b 2 ? ? 2 ( h1 x 1 ? h 2 x 2 ) ? 0
b1 ? b2

(2)对于线性方程组 ?

在*面上匹配两条直线 l1 和 l2,由题意

?1

?2

,x1(t0)>0,可将第一象限分为

三个区域。在最左边区域, x '1 , x ' 2 都大于 0;在中间区域, x '1 , x ' 2 都小于 0,在最右边
区域, x '1 , x ' 2 分别是大于 0 和小于 0.,由各区域中 x '1 , x ' 2 的取值可得到如下图形:

由图也可以看出,随着时间的增加,物种 1 最终能达到稳定值 终要灭亡。 4、蝴蝶效应与混沌解 考虑 Lorenz 模型
? x1 ( t ) ? ? ? x1 ( t ) ? x 2 ( t ) x 3 ( t ) ? ' x 2 (t ) ? ?? x 2 (t ) ? ? x3 (t ) ? ? x ' (t ) ? ? x (t ) x (t ) ? ? x (t ) ? x (t ) 1 2 2 3 ? 3
'

b1

? 1 h1

,物种 2 最

其中σ =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) 2(0)=0,x3(0) ,x ,重 复一的工作,看有什么现象发生。 解: (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 函数求解,可得数值解: Columns 1 through 9 0 0.1250 0.2500 0.3750 0.5000 0.5352 0.5705 0.6057 0.6409 0 0.0000 0.0000 -0.0000 -0.0000 0.0000 0.0000 0.0000 0.0000 0 0.0000 0.0000 -0.0000 -0.0000 0.0000 0.0000 0.0000 -0.0000

0.0000 -0.0000 0.0000 0.0000 Columns 10 through 18 0.6761 0.7114 0.9776 1.0105 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 Columns 19 through 27 1.0434 1.0763 1.2797 1.3186 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0001 0.0001 0.0000 0.0000 0.0002 0.0002 Columns 28 through 36 1.3575 1.3964 1.5656 1.5938 0.0000 0.0000 0.0000 0.0000 0.0002 0.0003 0.0021 0.0029 0.0004 0.0006 0.0045 0.0063

-0.0000

0.0000

0.0000

0.0000

-0.0000

0.7466 0.0000 0.0000 0.0000

0.7818 0.0000 0.0000 0.0000

0.8308 0.0000 0.0000 0.0000

0.8797 0.0000 0.0000 0.0000

0.9286 0.0000 0.0000 0.0000

1.1092 0.0000 0.0000 0.0000

1.1421 0.0000 0.0000 0.0000

1.1750 0.0000 0.0000 0.0000

1.2079 0.0000 0.0000 0.0001

1.2409 0.0000 0.0000 0.0001

1.4246 0.0000 0.0004 0.0008

1.4528 0.0000 0.0005 0.0012

1.4810 0.0000 0.0008 0.0016

1.5092 0.0000 0.0011 0.0023

1.5374 0.0000 0.0015 0.0032

............... Columns 5590 through 5598 99.5751 99.5921 99.6090 99.7069 99.7338 16.9457 16.5261 16.2010 17.1933 18.7894 -3.3551 -3.7119 -4.1098 -7.5438 -8.8677 -5.3476 -5.9274 -6.5941 -12.1944 -14.0725

99.6260 15.9854 -4.5568 -7.3519

99.6462 15.8978 -5.1636 -8.3766

99.6664 16.0348 -5.8601 -9.5370

99.6867 16.4476 -6.6527 -10.8209

Columns 5599 through 5607 99.7608 99.7878 99.8148 99.8940 99.9099 21.2186 24.4709 28.2629 36.7149 37.0528 -10.3044 -11.7459 -12.9932 -13.4302 -12.7919 -15.7415 -16.8309 -16.9274 -10.1013 -8.1005

99.8306 30.5354 -13.5361 -16.3648

99.8465 32.6533 -13.8800 -15.3302

99.8623 34.4645 -13.9854 -13.8707

99.8782 35.8466 -13.8351 -12.0805

Columns 5608 through 5613 99.9257 99.9416 99.9562 99.9708 99.9854 100.0000 36.8957 36.3239 35.5248 34.5435 33.4481 32.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) 2(0)不变,x3(0)=10-8 ,x 分别得到得到 x2~x1,x2~x3,x3~x1 的*面图如下:

可以看出,参数σ ,ρ ,β 值各减 1,初始值 x1(0) 2(0)不变,x3(0) ,x 数值变为 x3(0)=10-8,参数和初始值很小的改变,就会导致最后图形发生很大 的变化。

5、用微分方程考察共振现象 设物体沿 x 轴运动(如图所示)其*衡位置取为原点 0,物体的质量为 1,在 时间 t 物体的位置为 x(t)其所受的恢复为(如弹性力等)与物体所在位置的坐标成 正比,即 k2x,其中常数 k 称为恢复系数,运动过程所受的阻力(由于介质及摩擦 等)设与速度成正比,即 2 h
dx dt

,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)=B sinwt),建立相应的微分方程(该方程称为强迫振动方程).在上述参 数不变的情况下,取振幅 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)分别对上述图形讲行分析,并解释为什么会出现这些现象。
2

解: (1)根据 Newton 第二定律,F=ma,可得微分方程:m

d x dt
2

? m g ? k x ? 2h
2

dx dt

由题意知边界值:x(0)=1,x’(0)=0,令 y1=x,y2=
? d y1 ? y2 ? dt ? ? dy 可将二阶方程转化为: ? 2 ? m g ? k 2 y 1 ? 2 h y 2 ? dt ? y 1 ( 0 ) ? 1, y 2 ( 0 ) ? 0 ? ?

dx dt



其中,m=1;g=9.8;k=2 当 h=0 时,由 matlab 解得数值解: Columns 1 through 9 0 0.0000 0.0002 0.0002 1.0000 1.0000 1.0000 1.0000 Columns 10 through 18 0.0004 0.0006 0.0054 0.0108 1.0000 1.0000 1.0001 1.0003 Columns 19 through 27 0.0162 0.0216 0.2155 0.2956 1.0008 1.0014 1.1326 1.2461 ............ Columns 100 through 108 8.2321 8.3593 8.9852 9.0939 3.5019 3.2161 1.5214 1.3034 Columns 109 through 117 9.2025 9.3142 9.9124 10.0000 1.1393 1.0345 1.6412 1.8634 0.0000 1.0000 0.0000 1.0000 0.0000 1.0000 0.0001 1.0000 0.0001 1.0000

0.0009 1.0000

0.0011 1.0000

0.0022 1.0000

0.0032 1.0000

0.0043 1.0001

0.0271 1.0021

0.0541 1.0085

0.0812 1.0191

0.1083 1.0339

0.1353 1.0528

8.4615 2.9505

8.5636 2.6641

8.6658 2.3688

8.7679 2.0769

8.8766 1.7834

9.4260 1.0001

9.5377 1.0383

9.6494 1.1467

9.7371 1.2773

9.8247 1.4438

h=0 时,x(t)图形如下所示:

当 h=0.1 时,由 matlab 解得数值解: Columns 1 through 9 0 0.0000 0.0000 0.0002 0.0002 1.0000 1.0000 1.0000 1.0000 1.0000 Columns 10 through 18 0.0004 0.0006 0.0009 0.0054 0.0108 1.0000 1.0000 1.0000 1.0001 1.0003 Columns 19 through 27 0.0162 0.0216 0.0271 0.2156 0.2958 1.0008 1.0014 1.0021 1.1308 1.2417 ............. Columns 109 through 117 8.5631 8.6690 8.7749 9.3934 9.5289 2.5857 2.4563 2.3294 1.8880 1.8962

0.0000 1.0000

0.0000 1.0000

0.0001 1.0000

0.0001 1.0000

0.0011 1.0000

0.0022 1.0000

0.0032 1.0000

0.0043 1.0001

0.0541 1.0085

0.0812 1.0190

0.1083 1.0336

0.1353 1.0523

8.8809 2.2106

8.9868 2.1049

9.1223 1.9952

9.2579 1.9213

Columns 118 through 121 9.6467 9.7645 1.9353 2.0016

9.8822 2.0909

10.0000 2.1979

h=0.1 时,x(t)图形如下所示:

(2)如果还受到外力干扰,则微分方程变为:
m d x dt
2 2

? mg ? k x ? 2h
2

dx dt

? B s in w t

将其化为一阶方程组:
? d y1 ? y2 ? dt ? ? dy2 2 ? m g ? k y 1 ? 2 h y 2 ? B s in w t ? dt ? ? y 1 ( 0 ) ? 1, y 2 ( 0 ) ? 0 ? ?

其中,m=1;g=9.8;k=2;B=1;h=0.1; 用 Matlab 解得数值解: w=1 时 Columns 1 through 9 0 0.0000 0.0000 0.0000 0.0000 0.0001 0.0001

0.0002 0.0002 1.0000 1.0000 1.0000 1.0000 ............ Columns 118 through 121 9.8207 1.8973 9.8804 1.9210

1.0000

1.0000

1.0000

1.0000

1.0000

9.9402 1.9503

10.0000 1.9846

w=1.2 时 Columns 1 through 9 0 0.0000 0.0002 0.0002 1.0000 1.0000 1.0000 1.0000 ............ Columns 118 through 121 9.8281 1.6999 9.8854 1.7568 0.0000 1.0000 0.0000 1.0000 0.0000 1.0000 0.0001 1.0000 0.0001 1.0000

9.9427 1.8199

10.0000 1.8886

w=1.4 时 Columns 1 through 9 0 0.0000 0.0000 0.0000 0.0000 0.0002 0.0002 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 ............ Columns 118 through 125 9.6236 9.7309 9.8382 9.9454 9.9864 10.0000 2.2401 2.3171 2.4084 2.5103 2.5510 2.5647 w=1.4 时 Columns 1 through 9 0 0.0000 0.0000 0.0000 0.0000 0.0002 0.0002 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 ............ Columns 118 through 121 9.7823 9.8549 9.9274 10.0000

0.0001 1.0000

0.0001 1.0000

9.9591 2.5238

9.9727 2.5373

0.0001 1.0000

0.0001 1.0000

2.1150

2.0676

2.0285

1.9979

w=1.6 时 Columns 1 through 9 0 0.0000 0.0002 0.0002 1.0000 1.0000 1.0000 1.0000 ............ Columns 118 through 121 9.6933 9.7955 0.8062 0.7503

0.0000 1.0000

0.0000 1.0000

0.0000 1.0000

0.0001 1.0000

0.0001 1.0000

9.8978 0.7562

10.0000 0.8231

w=1.8 时 Columns 1 through 9 0 0.0000 0.0000 0.0000 0.0000 0.0002 0.0002 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 .............. Columns 118 through 121 9.7752 9.8501 9.9251 10.0000 0.8629 1.0850 1.3383 1.6174

0.0001 1.0000

0.0001 1.0000

w=2.0 时 Columns 1 through 9 0 0.0000 0.0002 0.0002 1.0000 1.0000 1.0000 1.0000 ............. Columns 118 through 121 9.7538 9.8359 2.3278 2.6023 ............. w=3.0 时 Columns 1 through 9 0 0.0000 0.0002 0.0002 1.0000 1.0000 1.0000 1.0000 ...........

0.0000 1.0000

0.0000 1.0000

0.0000 1.0000

0.0001 1.0000

0.0001 1.0000

9.9179 2.8705

10.0000 3.1243

0.0000 1.0000

0.0000 1.0000

0.0000 1.0000

0.0001 1.0000

0.0001 1.0000

Columns 118 through 121 9.7495 9.8330 2.2414 2.3303

9.9165 2.4145

10.0000 2.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.9 5.3 7.2 9.6 12.9 17.1 23.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'); hold on 在以上每一次求时顺便画出图形,得到最后的总图如下

(2)利用 Logistics 模型,将以上所有数据进行拟合,编写程序如下 t = 0 : 10 : 190; x = [3.9 5.3 7.2 9.6 12.9 17.1 23.2 31.4 38.6 50.2 62.9 76.0 92.0 106.5 123.2 131.7 150.7 179.3 204.0 226.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 所示),建立反映肥料生成机理的数学 模型,提出改善肥料生产的建议。

解答: 1.问题条件假设 1)每次堆肥的质量相等 2)所给的几次堆肥混合物的比例仅由当天的实际情况决定 3)所有分离肥仓工作条件相同 4)每磅蔬菜可提供的氧气量相等 5)细菌消耗的溶解氧气完全由蔬菜叶提供 6)每天提供的废物混合物中的化学成分大致相同 7)废物混合物在喂给细菌前混合均匀并保持良好的通气环境。 2.问题分析 堆肥是利用微生物的分解作用将有机废物转化成无害稳定形式的生物化学过程, 要提高堆肥产量方法之一就是通过增强细菌的生长繁殖能力以提高分解率。 细菌群体的增长一般要经历延滞期,加速生长期、对数生长期、减数生长期、稳 定期、加速死亡期和对数死亡期。另外,对数生长期培养基中所有养分都过剩, 细菌可以充分繁殖, 其倍增速率恒定, 取决于底物浓度、 温度、 水活度、 供氧量。 对于当前该餐厅来说底物浓度由每天的剩余食物、蔬菜叶和碎纸屑决定。碎纸屑 是吸收水分的调理剂, 微菌呼吸所需要的溶解氧由蔬菜也提供,水活度可以通过

测定相对湿度来决定,其关系式是: 相对湿度:B=P/P0*100% 水活度: aw=P/P0 其中 P 为该溶液蒸汽压,P0 为纯水蒸汽压, 在这个堆肥系统中,可供微菌消耗的营养底物和溶解氧都是有限的,他们的消耗 会对微菌生长率产生重要影响, 一种微菌消耗营养底物的速率和底物浓度之间的 关系式为: ds/dt=—KmSX/(Ks+S) ,式中 ds/dt 表示为底物的有效消耗率,x 表示为微菌浓 度,Km 表示为最大有效系数,在高浓度营养底物中最大的物料消耗率(物质质 量/微菌一天的质量) ;Ks 表示为半速系数(质量/体积) 表示为有限底物的浓 ,S 度(质量/体积) 微菌生长过程是一个生物化学过程,其生长率和温度之间满足公式:K=Ae^(— Ea/RT),式中 K 表示为反应速度常数,T 表示为反应的绝对温度,R 表示气体常 数,Ea 表示反应活化能,A 表示频率分子。 我们知道, 各种微菌都有一个最适生长温度,如果温度控制在最适值是微菌生长 速率最高。 微菌生长对水活度也有一定的要求,与微菌最高生长速率相对应的有 一最适水活度。 优化堆肥意味着尽可能少的时间内生产出高质量的肥料,参数的 优化依赖于所应用的系统。 3 模型的建立 假设每天投入的废料比是随意的,仅仅取决于当天的情况,首先从已知数据中得 到经验最佳比例, 由于假设各次堆肥后肥料质量相同,因此堆肥时间较短就对应 了较好的废料配比,将 12 组数据按其堆肥日期及完全堆肥时间分成三组,每组 中的较优比例见下表: 分组 每组中的最短天数 比例 1~4 No.4: 26 203:82:0 5~8 No.5: 33 79:28:0 9~12 No.9: 49 82:44:9 上述模型显然过于粗糙 模型二: 营养底物和氧气是细菌生长的两种底物,物耗公式为: dSi/dt=—KmiXSi/(Ksi+Si) (1) i=1 时代表营养底物中可降解得有机物。i=2 时代表氧气。 1)在高浓度底物中,物料的转化过程很迅速,进一步增加底物浓度就不再引起 底物转化率的提高,可假设 S1>>Ks1,S2>>Ks2,则(1)式简化为:dSi/dt*(1/X)= —Kmi, (2) 这是关于底物浓度的零级反应,消去 X 得到:dS1=Km1/Km2*dS2 (3) 积分后得:S1(t)=Km1/Km2*S2(t)+{S1(0)—Km1/Km2*S2(0)} (4) 2)在低浓度底物中,可假设 S1<<Ks1,S2<<Ks2,则(1)式可简化为: dSi/dt*(1/X)=—Kmi/Ksi*Si (5) 消去 X 得:dS1/S1 *1/X=Km1*Ks2/(Km2*Ks1)* dS2/S2 (6) 积分后得:S1(t)=S1(0)[S2(t)/s2(0)]^[Km1*Ks2/(Km2*Ks1)] (7) 3)当 Si=Ksi,i=1,2 时有:dS1=Km1/Km2*dS2 下面的步骤就和 1 相同了 4 结果分析: 由上面的计算可知,当 S>>Ks 时,简化方程为(4)且 Km1=Ae^(—Ea/RT),其 中,

R=8.31J/(mol*K) ,Ea=51.3KJ/mol T=25°C=298K 时,Km=12 克底物/(克微菌*天) 可确定:A=12e^(51.3*1000/8.31*298) 当 T=49°C=322K 时,可得:Km1=A*e^(51.3*1000/8.31*322)=55.7 不 同 微 生 物 的 耗 氧 速 度 不 同 , 这 里 可 取 86 , 经 过 单 位 换 算 , 可 估 算 出 Km2=86*0.001*24*223/22.4=20.548 从而 Km1/Km2=55.7/20.548=2.711 下面从已知数据中选取几组典型的数来定性说明(4)式在不同情况下的物理意 义, 设废物浆中含有 40%的有机挥发性物质,有机物的降解为 66%,并设蔬菜叶中 含有 25%的氧气, A)对于第二组数据:S1(0)=112*40%*66%=22.704 S2(0)=79*25%=19.75, S1(t)=2.711S2(t)—23.974,当 S1(t)=0 时,有 S2(t)=8.84,这表明在此情形下废物 浆中所含有的可降解的有机物可以全部被降解,氧气的量是充足的。 B)对第五组数据:S1(0)=20.856 S2(0)=7,S1(t)=2.711S2(t)+1.879,在时间足 够的情况下,S2(t)=0,S1(t)=1.879,这表明氧气的量是不充足的,废物浆中 所含可降解的有机物最多只能有 S1(0)-S1(∞)/S1(0)=91%被降解 由此我们可以确定出废物浆和蔬菜叶的最优比例关系, 最优情况是在某个时刻 t0 有 S1(t0)=S2(t0)=0,即废物浆中可降解的有机物和蔬菜叶所提供的氧气据 别利用完。设堆肥原料中废物浆的量为 x,蔬菜叶的量为 y,则由 S1(t0)=S2 (t0) 可得: =0 S1(0)—2.711S2(0)=0,故有 x*40%*66%=2.711*y*25%,从而得到 x: y=2.567:1 对 S<<Ks 及 S=Ks 的两种情况,可用类似的方法进行分析。 5 模型评价 我们将堆肥过程中的微菌与堆肥原料之间复杂的相互作用, 根据三种不同的情形 进行了较为合理的简化, 建立了一个比较简单的反映废物浆和蔬菜叶关系的数学 模型, 并对其中两种情况做出了物理解释,但没有考虑微菌与堆肥原料之间更为 复杂的相互促进和限制作用,模型建立过于简单,另外,由于所给数据有限及堆 肥问题的复杂性, 对一些影响堆肥的因素只做了定性分析,而缺乏更精确地定量 分析, 6 实施建议 1)堆肥过程中,每隔一段时间添加一定量的新鲜废物,并从堆肥仓中排除等量 的液体,使营养底物的浓度保持在一定水*之上,促进微菌大量繁殖。 2)堆肥过程中要保证适量的水分,含水量保持在原料湿度 60%-70%为宜 3) 堆肥初期应经常通气以利于微菌活动,后期则应少通气以利于腐殖后形成和 减少养分损失 4)微菌为了进行呼吸和繁殖,对堆肥的碳氮比有一定要求,多数微菌碳氮比为 25:1 5)控制堆温使微菌能长时间在最适温度下生存 6)中性和微碱性条件有利于微菌活动,可加入一些石灰、草木灰等碱性物质 7)具体操作时,应对堆肥进行充分搅拌以使混合均匀




友情链接: