数学实验2

数学实验.实验2

4.某日凌晨一住所发生一件凶杀案，警方于6时到达现场后测得尸温26℃，室温17℃，2小时后尸温下降了3℃，试根据冷却定理建立差分方程，估计凶杀案发生的时间（可设正常体温为37℃）

解：设t0为正常体温，y为死亡后人的体温，x为时间。差分方程由冷却定理得。

冷却定律：

T=C+(T0-C)*e^-kt

C是室温，T0是人体常温37℃，T为死亡后第t小时的体温，k为可求常数。

程序：

函数

function y=swsj1(x)

%带参数求得k=0.2027，令a=k

t0=37;

a=0.2027

y=17+(t0-17)*exp(-a*x)

end

n=6;

x=0:n;

y=swsj1(x)

plot(x,y','r*')

gtext('死亡时间'),

运行结果：

a =

0.2027

y =

37.0000 33.3305 y =

30.3342 27.8877 25.00 24.25 22.9271

37.0000 33.3305 30.3342 27.8877 25.00 24.25 22.9271

由结果知此人的死亡时间为2 :00.

5.据报道，某种山猫在较好、中等及较差的自然环境下，年平均增长率分别是1.68%，0.55%和-4.50%，假定开始时有100只山猫，按以下情况讨论山猫数量逐年变化过程及趋势：

（1）3种自然环境下25年的变化过程（作图）；

（2）如果每年捕获三只，会发生什么情况？山猫会灭绝吗？如果每年捕获一只呢？

（3）在较差的自然环境下，如果想使山猫数量稳定在60只左右，每年需要人工繁殖多少只？

（1）3种自然环境下25年的变化过程（作图）

程序:

function x=cat(x0,n,r)

x=x0;

for k=1:n

x(k+1)=(1+r)*x(k);

end

k=(0:25)';

y1=cat(100,25,0.0168);

y2=cat(100,25,0.0055);

y3=cat(100,25,-0.0450);

round([k,y1',y2',y3'])

plot(k,y1,k,y2,'r:',k,y3,'g--'),grid

gtext('较好环境'),gtext('中等环境'),gtext('较差环境')

结果：

ans =

0 100 1 102 2 103 3 105 4 107 5 109 6 111 7 112 8 114 9 116 100

96

91

87

83

79

76

72

69

66

100 101 101 102 102 103 103 104 104 105

10 118 106 63

11 120 106 60

12 122 107 58

13 124 14 126 15 128 16 131 17 133 18 135 19 137 20 140 21 142 22 144 55

52

50

48

46

44

42

40

38

36

107 108 109 109 110 110 111 112 112 113

23 147 113 35

24 149 114 33

25 152 115 32

（2）如果每年捕获三只，会发生什么情况？山猫会灭绝吗？如果每年捕获一只呢？

（i）每年捕获三只

程序:

function x=cat(x0,n,r)

x=x0;

for k=1:n

x(k+1)=(1+r)*x(k)-3;

end

command窗口程序同第一问

结果：

ans =

0 100 100 100

1 99 98 93

2 97 95 85

3 96 93 78

4 95 90 72

5 93 88 66

6 92 85 60

7 90 83

8 80 49

9 87 77 43

10 86 11 84 12 83 13 81 14 79 15 78 16 76 17 74 18 73 19 71 75 39

72 34

70 29

67 25

21

62 17

59 13

56 10

6

51 3

20 69 48 0

21 67 46 -3

22 65 43 -6

23 63 40 -9

24 61 37 -11

25 59 35 -14

(ii) 每年捕获一只

程序:

function x=cat(x0,n,r)

x=x0;

for k=1:n

x(k+1)=(1+r)*x(k)-1;

end

command窗口程序同第一问

结果:

ans =

0 100 100 100

1 101 100 95

2 101 99

3 102 99 84

4 103 98 79

5 104 98 75

6 104 97 70

7 105 97 66

8 106 9 107 10 107 11 108 12 109 13 110 14 111 15 111 16 112 17 113 96 62

96 59

95 55

95 51

94 48

94 45

93 42

93 39

92 36

92 34

18 114 92 31

19 115 91 29

20 116 91 26

21 117 22 118 23 119 24 120 25 121 90 24

90 22

20

88 18

88 16

答:若每年捕获三只,无论在什么环境情况下,均会灭绝;若每年捕获一只,则在中等和较差的环境下会灭绝.

（3）在较差的自然环境下，如果想使山猫数量稳定在60只左右，每年需要人工繁殖多少只？

程序:

syms b

solve(60-(1-0.0450)*60-b)

结果:

ans =

27/10

答: 在较差的自然环境下，如果想使山猫数量稳定在60只左右，每年需要人工繁殖3只.

8.在某种环境下猫头鹰的主要食物来源是田鼠，设田鼠的年平均增长率为r1，猫头鹰的存在引起的田鼠增长率的减少与猫头鹰数量成正比，比例系数a1；猫头鹰的年平均减少率为r2；田鼠的存在引起的猫头鹰减少率增加与田鼠的数量成正比，比例系数为a2，建立差分方程模型描述田鼠和猫头鹰共处时数量变化规律，对以下情况作图给出50年的变化过程。

解：设x为田鼠数，y为猫头鹰数。

程序如下：

函数：

function z=ms(x0,y0,r1,r2,a1,a2,n)

x=x0;

y=y0;

for k=1:n

x(k+1)=(1+r1-a1*y(k))*x(k);

y(k+1)=(1-r2+a2*x(k))*y(k);

end

z=[x',y'];

运行结果；

ans =

37.2001

z =

2

y =

0

37.2001

33.5035

30.4834

28.0159

26.0000

（1）设r1=0.2，r2=0.3,a1=0.001，a2=0.002,开始时有100只田鼠和50只猫头鹰。

z=ms(100,50,0.2,0.3,0.001,0.002,50);

k=(0:50)';

[k,z]

plot(k,z(:,1),k,z(:,2)),grid,

gtext('x(k)'),gtext('y(k)'),

(2)r1,r2,a1,a2同上，开始时有100只田鼠和200只猫头鹰。

z=ms(100,200,0.2,0.3,0.001,0.002,50);

k=(0:50)';

[k,z]

plot(k,z(:,1),k,z(:,2)),grid,

gtext('x(k)'),gtext('y(k)'),

(3)适当改变参数a1,a2（初始值同上）。

令a1=0.003,a2=0.001,有如下数据图像

z=ms(100,200,0.2,0.3,0.003,0.001,50);

k=(0:50)';

[k,z]

plot(k,z(:,1),k,z(:,2)),grid,

gtext('x(k)'),gtext('y(k)'),

（4）求差分方程的平衡点，它们稳定吗？

令x(k)=x(k+1)=x,y(k)=y(k+1)=y,得平衡点是(150，200)。

但是由图像可以看出它们并不稳定。

9.研究将鹿群放入草场后草和鹿两种群的相互作用，草的生长遵从Logistic规律，年固有增长率0.8，最大密度为3000（密度单位），在草最茂盛时每只鹿每年可吃掉1.6（密度单位）的草，若没有草，鹿群的年死亡率高达0.9，而草的存在可使鹿的死亡的以补偿，在草最茂盛时补偿率为1.5，作出一些简化假设，用差分方程模型描述草和鹿两种群数量的变化过程，就一下情况进行讨论，

（1）比较将100只鹿放入密度为1000和密度3000的草场两种情况。

（2）适当改变参数，观察变化趋势。

解：设x为鹿群数，y为草的数量，H为草场的最大密度。

（1）

程序：

function z=chl(x0,y0,n)

H=3000;

x=zeros(1,n);y=x;

a=x;b=x;

a(1)=1.6;b(1)=1.5;

x(1)=x0;

y(1)=y0;

for k=1:n

y(k+1)=y(k)+0.8*(1-y(k)/H)*y(k)-a(k)*x(k);

x(k+1)=x(k)+(b(k)-0.9)*x(k);

a(k+1)=a(1)*(y(k+1)/H);

b(k+1)=b(1)*(y(k+1)/H);

end

z=[x',y'];

end

z1=chl(100,3000,60)

z2=chl(100,1000,60)

k=0:60;

plot(k,z1(:,1),'b',k,z1(:,2),'b',k,z2(:,1),'r',k,z2(:,2),'r'),grid,

gtext('x，“鹿群”'),gtext('y“草”'),

运行结果：

z1 =

1.0e+003 *

0.1000 3.0000

0.1600 2.8400

0.2432 2.7188

0.39 2.5700

0.4916 2.3782

0.6337 2.1490

0.7443 1.9104

0.78 1.7071

0.74 1.5806

0.6668 1.76

0.5826 0.5234 0.4967 0.5001 0.5268 0.5674 0.6100 0.23 0.6558 0.91 1.5966

1.6980

1.8136

1.9070

1.92

1.9501

1.9059

1.8419

1.7797

1.73

0.6285 1.7203

0.6035 1.7307

0.5826 1.7595

0.5708 0.5693 0.5765 0.58 0.6021 0.6122 0.6169 0.6157 0.6100 0.6026 1.7949

1.8253

1.8429

1.8449

1.8338

1.8152

1.7960

1.7817

1.7755

1.7776

0.5958 1.7858

0.5916 1.7965

0.5906 1.8062

0.5924 0.5961 0.6001 0.6034 0.6051 0.6049 0.6034 0.6012 0.5990 0.5976 1.8123

1.8137

1.8109

1.8055

1.7996

1.7949

1.7926

1.7929

1.7952

1.7984

0.5971 1.8015

0.5976 1.8035

0.5986 1.8042

0.5999 0.6009 0.6015 0.6015 0.6011 0.6005 0.5998 0.5993 0.5991 0.5992 1.8035

1.8019

1.8001

1.7986

1.7978

1.7978

1.7984

1.7994

1.8003

1.8010

0.5995 1.8013

0.5999 1.8011

z2 =

1.0e+003 *

0.1000 1.0000

0.1600 1.3733

0.1259 1.8519

0.1291 2.2945

0.1611 2.5682

0.2229 2.33

0.3169 2.5805

0.4406 2.4330

0.5800 2.2292

0.7045 1.9978

0.7742 1.7811

0.7668 1.6246

0.6996 0.6143 0.50 0.5061 0.4983 0.5163 0.5518 0.5940 0.6307 0.6516 1.5560

1.5746

1.6573

1.76

1.8721

1.9377

1.9531

1.9236

1.8663

1.8027

0.6525 1.7518

0.6368 1.7252

0.6130 1.7258

0.5902 0.5749 0.5695 0.5735 0.5841 0.5971 0.6086 0.61 0.61 0.6123 1.7480

1.7813

1.8141

1.8368

1.8447

1.8384

1.8224

1.8032

1.7868

1.7775

0.60 1.7765

0.5983 1.7825

0.5931 1.7924

0.5908 0.5916 0.5946 0.5986 0.6022 0.6045 0.6051 0.6040 0.6020 0.5998 1.8027

1.8102

1.8134

1.8121

1.8076

1.8018

1.7966

1.7934

1.7927

1.7942

0.5981 1.7972

0.5973 1.8003

0.5974 1.8028

0.5982 0.5994 0.6005 0.6013 0.6015 0.6013 0.6007 0.6000 0.5995 0.5992 1.8040

1.8038

1.8026

1.8008

1.7992

1.7981

1.7978

1.7981

1.7990

1.8000

由图可知，鹿群与草的数量符合寄生模型，最终鹿和草的数量都趋于一个稳定的值，由差分方程（y(k+1)=y(k)+0.8*(1-y(k)/H)*y(k)- a(1)*(y(k+1)/H)*x(k);

x(k+1)=x(k)+(b(1)*(y(k+1)/H)-0.9)*x(k);

令x=x（i）,y=y(i),i=1,2…n)

可算出稳定值约为(600,1800)。

（2）

function z=chl(x0,y0,a1,b1,c,d,n)

H=3000;

x=zeros(1,n);y=x;

a=x;b=x;

a(1)=a1;b(1)=b1;

x(1)=x0;

y(1)=y0;

for k=1:n

y(k+1)=y(k)+c*(1-y(k)/H)*y(k)-a(k)*x(k);

x(k+1)=x(k)+(b(k)-d)*x(k);

a(k+1)=a(1)*(y(k+1)/H);

b(k+1)=b(1)*(y(k+1)/H);

end

z=[x',y'];

end

z1=chl1(100,3000,1.6,1.8,0.8,0.9,60)

z2=chl1(100,1000,1.6,1.8,0.8,0.9,60)

k=0:60;

plot(k,z1(:,1),'b',k,z1(:,2),'b',k,z2(:,1),'r',k,z2(:,2),'r'),grid,

gtext('x“鹿群”'),gtext('y“草”'),

运行结果：

z1 =

1.0e+003 *

0.1000 3.0000

0.1900 2.8400

0.3428 2.6734

0.5841 2.4175

0.9056 2.0400

1.1990 1.5769

1.23 1.1669

1.0037 0.9567

0.6765 0.4597 0.3608 0.3450 0.3970 0.5190 0.7108 0.9365 1.0968 1.0783 0.9659

1.1413

1.4272

1.7512

2.0121

2.1162

2.0292

1.7852

1.4718

1.2106

0.11 1.0921

0.6730 1.1287

0.5231 1.2868

0.4562 0.4605 0.5287 0.6561 0.8224 0.9715 1.0224 0.9383 0.7777 0.6308 1.5157

1.7469

1.9016

1.9224

1.8022

1.5874

1.3629

1.2147

1.1851

1.2671

0.27 1.42

0.5187 1.6121

0.5536 1.7628

0.09 0.76 0.06 0.9596 0.9313 0.8238 0.6994 0.6083 0.5688 0.5815 1.8239

1.7725

1.6291

1.4509

1.3077

1.2483

1.2829

1.3918

1.5371

1.6704

0.10 1.7446

0.7350 1.7323

0.8375 1.6388

0.9072 0.9081 0.8401 0.7418 0.6572 0.6108 0.6086 0.78 0.7187 0.8017 1.5017

1.3751

1.3049

1.3101

1.3822

1.4940

1.6073

1.6825

1.6924

1.6338

0.8660 1.5305

0.8819 1.4233

0.8413 1.3523

0.7667 1.3398

z2 =

1.0e+003 *

0.1000 1.0000

0.1900 1.3733

0.1756 1.8299

0.2103 2.2295

0.3024 2.4375

0.4725 2.4101

0.7304 2.1819

1.0293 1.8079

1.2195 1.3902

1.1391 1.0828

0.80 0.5868 0.4320 0.3742 0.38 0.4728 0.6248 0.8294 1.0209 1.01 0.9786

1.0604

1.2770

1.5695

1.8550

2.0357

2.0458

1.8846

1.6115

1.3308

0.9786 1.1501

0.7731 1.1172

0.5956 1.2175

0.4946 0.4678 0.5058 0.6037 0.7497 0.9058 1.0017 0.9766 0.8469 0.6943 1.4095

1.6355

1.8226

1.9032

1.8471

1.67

1.4582

1.2787

1.1996

1.2337

0.5834 1.3580

0.5336 1.5301

0.33 1.6944

0.6066 0.7134 0.8376 0.9314 0.9447 0.8692 0.7510 0.71 0.5876 0.5790 1.7934

1.7902

1.6866

1.5239

1.3667

1.2734

1.2694

1.3468

1.4757

1.6131

0.6183 1.7116

0.6968 1.7352

0.7951 1.6756

0.87 0.90 0.8686 0.7818 0.6919 0.6308 0.6112 0.6332 0.6904 0.7679 1.5568

1.4262

1.3334

1.3083

1.3530

1.4480

1.5601

1.6506

1.6871

1.6566

0.8401 1.5716

0.8761 1.4661

z1=chl1(100,3000,1,1.5,0.8,0.85,60)

z2=chl1(100,1000,1,1.5,0.8,0.85,60)

k=0:60;

plot(k,z1(:,1),'b',k,z1(:,2),'b',k,z2(:,1),'r',k,z2(:,2),'r'),grid,

gtext('x“鹿群”'),gtext('y“草”'),

改变参数可看出a1明显影响着两者数量的平衡性；b1明显影响其震荡的频率及次数。