MATLAB求差分方程模型

x1 (k 1) 0.6 x1 (k ) 0.2 x2 (k ) 0.1x3 ( k ) x2 (k 1) 0.3 x1 ( k ) 0.7 x2 ( k ) 0.3 x3 ( k ) x (k 1) 0.1x ( k ) 0.1x ( k ) 0.6 x ( k ) 1 2 3 3
c2
k 2
1,2 1, xk 0(k )
1,2 1, xk ( k )
• 本例中，用待定系数的方法可以求出 b=0.18时,c1=95.64, c2=4.36 ， (1, 2 ) (0.9430, 0.0430) k k 这样 x 95.64(0.9430) 4.36(0.0430)
k
5 实际上， 1,2
10 2
b
插入的数 的大小， 其必须满 足一个递 减趋势， 增长值为 负
植物能一直繁殖下去的条件是b>0.191
线性常系数差分方程组
• 汽车租赁公司的运营
一家汽车租赁公司在3个相邻的城市运营，为方 便顾客起见公司承诺，在一个城市租赁的汽车可以在 任意一个城市归还。根据经验估计和市场调查，一个 租赁期内在A市租赁的汽车在A,B,C市归还的比例分别 为0.6,0.3,0.1;在B市租赁的汽车归还比例0.2,0.7,0.1;C 市租赁的归还比例分别为0.1,0.3,0.6。若公司开业时将 600辆汽车平均分配到3个城市，建立运营过程中汽车 数量在3个城市间转移的模型，并讨论时间充分长以 后的变化趋势。
0.1 A 0.6 0.7 0.3 B C
A
0.2
B
0.1
C
0.6 C
假设在 每个租 赁期开 始能把 汽车都 租出去， 并都在 租赁期 末归还
A
B 0.1
0.3
模型及其求解
• 记第k个租赁期末公司在ABC市的汽车数量 分别为x1(k),x2(k),x3(k)（也是第k+1个租赁 期开始各个城市租出去的汽车数量），很 容易写出第k+1个租赁期末公司在ABC市的 汽车数量为(k=0,1,2,3· · · )
Xk= a1bcXk-1 + a2b(1-a1)bcXk-2
• 实际上，就是Xk= pXk-1 + qXk-2 我们需 要知道x0,a1,a2,c, 考察b不同时，种子繁 殖的情况。在这里假设 • X0=100,a1=0.5,a2=0.25,c=10,b=0.18~0.20 • 这样可以用matlab计算了
当a=3， a=4时，如图
100 99.5 99 98.5 98 97.5 97 96.5 96
0
2
4
6
8
10
12
14
16
18
20
112
110
108
106
104
102
100
0
2
4
6
8
10
12
14
16
18
20
高阶线性常系数差分方程 ——一年生植物的繁殖
• 一年生植物春季发芽，夏天开花，秋季产种， 没有腐烂，风干，被人为掠取的那些种子可以 活过冬天，其中一部分能在第2年春季发芽， 然后开花，产种，其中的另一部分虽未能发芽， 但如又能活过一个冬天，则其中一部分可在第 三年春季发芽，然后开花，产种，如此继续， 一年生植物只能活1年，而近似的认为，种子 最多可以活过两个冬天，试建立数学模型研究 这种植物数量变化的规律，及它能一直繁殖下 去的条件。
作图观察5年以后数量的变化趋势：在指令窗口输入 >> y1=qczl(5); 300 >> k=0:5; 280 >> plot(k,y1)
260 240 220 200 180 160 140 120
0
0.5
1
1.5
2
2.5
3
3.5
4
4.5
5
练习： 1. 若最初600辆汽车平均分配在三个城市， 作图观察30年以后数量的变化趋势 2. 若最开始600辆汽车都在A市，变化时间 充分长以后，各城市汽车数量的变化情况 如何？
1. 在指令窗口输入： >>X0=100; >> r=0.0194; >> n=20; >>Xn=(1+r)^n*X0 • 执行后得到： Xn =146.8563
X0
2. 建立M脚本文件 在M文件编辑器中输入： X0=100; r=0.0194; n=20; Xn=(1+r)^n*X0 单击工具栏中的保存按钮，文件名sqh1.m 保存以后选择Debug:run菜单，则会在指令窗 口输出： >> Xn =146.8563
function Xn=sqh(n,r,X0) Xn=(1+r)^n*X0; end 单击工具栏中的保存按钮，文件名sqh3.m 保存以后在指令窗口调用sqh3函数,输入： >>Xn＝sqh3(20,0.0194,100) 输出同上。
Matlab实现
Xk+1 ＝(1+r)Xk
function x=sqh(n,r) x(1)=100; for k=1:n x(k+1)=(1+r)*x(k); end 单击工具栏中的保存按钮，文件名sqh4.m 保存以后在指令窗口调用sqh4函数,输入： >>xn＝sqh4(20,0.0194) 输出同上。
设在中等自然条件下每年孵化5只鹤 放入保护区. >> k=0:20; %一个行向量 >> y=rgsqh(20,-0.0324,5); >> plot(k,y)
130
125
120
115
110
105
100
0
2
4
6
8
10
12
14
16
18
20
每年放养多少只会保证物种数量的稳定？
令Xk=Xk+1 得a=－rX0 当X0＝100，r=－0.0324时， a=3.24
0
2
4
6
8
10
12
14
16
18
20
在同一坐标系下画图
>> k=0:20; %一个行向量 >> y1= sqh4(20, 0.0194); >> y2= sqh4(20, -0.0324); >> y3= sqh4(20, -0.0382); >> plot(k,y1,k,y2,':',k,y3,'r')
• 用矩阵表示
x1 (k 1) 0.6 0.2 0.1 x1 ( k ) x ( k 1) 0.3 0.7 0.3 x ( k ) 2 2 x (k 1) 0.1 0.1 0.6 x (k ) 3 3
3.建立M函数文件 格式： function 输出变量＝函数名称（输入变量） 例如：输入 function Xn=sqh(n,r) X0=100; Xn=(1+r)^n*X0; end 单击工具栏中的保存按钮，文件名sqh2.m 保存以后在指令窗口调用sqh2函数,输入： >>Xn＝sqh2(20,0.0194) 输出同上。
160
140
120
100
80
60
40
0
2
4
6
8
ቤተ መጻሕፍቲ ባይዱ
10
12
14
16
18
20
• 人工孵化是挽救濒危物种的措施之一，如果 每年孵化a只鹤放入保护区，观察在自然条件 下沙丘鹤的数量如何变化
Xk+1=（1+r）Xk +a
在M文件中输入： function X=rgsqh(n,r,a) X(1)=100; for k=1:n X(k+1)=（1＋r)*X(k)+a; end
记一棵植物春季产种的平均数为c,种子能活过一个冬 天的(1岁种子)比例为b,活过一个冬天没有发芽又活过一 个冬天的（2岁种子）比例仍为b,1岁种子发芽率a1，2岁 种子发芽率a2。
• 设c,a1,a2固定,b是变量,考察能一直繁殖的条件
• 记第k年植物数量为Xk,显然Xk与Xk-1,Xk-2有 关，由Xk-1决定的部分是 a1bcXk-1,由Xk-2决 定的部分是 a2b(1-a1)bcXk-2 Xk= a1bcXk-1 + a2b(1-a1)bcXk-2
结果分析：Xk＋pXk-1 + qXk-2＝0 (1)
x1+px0=0
k
(2)
• 对高阶差分方程可以寻求形如 xk 2 p q 0 的解。代入(1)式得 称为差分方程的特征方程。差分方程的特征根：
1,2
p
p 2 4q 2
k 1 1
方程(1)的解可以表为 xk c C1,C2 由初始条件x0,x1确定。
观察n年以后的3个城市的汽车数量变化情况
1. 直接在命令窗口（commond window）作
>> A=[0.6,0.2,0.1;0.3,0.7,0.3;0.1,0.1,0.6]; >> n=10; >> for k=1:n x(:,1)=[200,200,200]'; x(:,k+1)=A*x(:,k); end >> x(:,k+1)
输出： ans = 179.9324 299.9895 120.0781
2. 建立M文件 function x=qczl(n) A=[0.6,0.2,0.1;0.3,0.7,0.3;0.1,0.1,0.6]; x(:,1)=[200,200,200]'; for k=1:n x(:,k+1)=A*x(:,k); end
用Matlab求解差分方程问题
一阶线性常系数差分方程
高阶线性常系数差分方程
线性常系数差分方程组
一阶线性常系数差分方程
Florida沙丘鹤属于濒危物种，它在较好自然环境下， 年均增长率仅为1.94%，而在中等和较差环境下年均增 长率分别为 -3.24% 和 -3.82%，如果在某自然保护区内 开始有100只鹤，建立描述其数量变化规律的模，并作 数值计算。 模型建立: 记第k年沙丘鹤的数量为Xk,年均增长率为r， 则第k+1年鹤的数量为
Xk= a1bcXk-1 + a2b(1-a1)bcXk-2
假设 X0=100,a1=0.5,a2=0.25,c=10,b=0.18~0.20 简化为 Xk= pXk-1 + qXk-2 • • • • • • • • function x=zwfz(x0,n,b) c=10;a1=0.5;a2=0.25; p=a1*b*c;q=a2*b*(1-a1)*b*c; x(1)=x0; x(2)=p*x(1); for k=3:n x(k)=p*x(k-1)+q*x(k-2); end
>> k=0:20; >> y1=zwfz(100,21,0.18); >> y2=zwfz(100,21,0.19); >> y3=zwfz(100,21,0.20); >> plot(k,y1,k,y2,':',k,y3,'r')
250
200
150
100
50
0
0
2
4
6
8
10
12
14
16
18
20
141.3199 144.0615 146.8563
利用plot 绘图观察数量变化趋势
>>k=0:20； >>y1= sqh4(20,0.0194); >>plot(k,y1) 得到Figure1 在Edit菜单下Copy Figure，可以将图 粘贴到Word文档。
150 145 140 135 130 125 120 115 110 105 100
xn = Columns 1 through 6 100.0000 101.9400 103.9176 105.9336 107.9888 110.0837 Columns 7 through 12 112.2194 114.3964 116.6157 118.8780 121.1843 123.5353 Columns 13 through 18 125.9318 128.3749 130.8654 133.4042 135.9922 138.6305 Columns 19 through 21
35
得
Xk+1 ＝(1+r)Xk , k=0,1,2· · · · · · Xn ＝ (1+r)^n X0
已知X0=100, 在较好，中等和较差的自然环境下 r=0.0194 , -0.0324 和 -0.0382 我们利用Matlab编程， 递推20年后观察沙丘鹤的数量变化情况
Matlab实现 Xn ＝ (1+r)^n
1. 可以看到时间充分长以后3个城市汽车 数量趋于180，300，120 2. 若最开始600辆汽车都在A市，可以看 到变化时间充分长以后，各城市汽车数 量趋于稳定，与初始值无关
利用线性规划知识；找到凸集，因为所谓的最优解，一般都在凸集
600
500
400
300
200
100
0
0
5
10
15
20
25
30