用matlab求解差分方程

用Matlab求解差分方程问题一阶线性常系数差分方程高阶线性常系数差分方程线性常系数差分方程组差分方程是在离散时段上描述现实世界中变化过程的数学模型•例1、某种货币1年期存款的年利率是r，现存入M元，问年后的本金与利息之和是多少？•Xk+1=(1+r)xk,k=0,1,2·····•以k=0时x0=M代入，递推n次可得n年后本息为1nnxrM•污水处理厂每天可将处理池的污水浓度降低一个固定比例q，问多长时间才能将污水浓度降低一半？•记第k天的污水浓度为ck,则第k+1天的污水浓度为ck+1=(1-q)ck,k=0,1,2,····从k=0开始递推n次得以cn=c0/2代入即求解。0(1)nncqc一阶线性常系数差分方程•濒危物种的自然演变和人工孵化•问题Florida沙丘鹤属于濒危物种，它在较好自然环境下，年均增长率仅为1.94%，而在中等和较差环境下年均增长率分别为-3.24%和-3.82%，如果在某自然保护区内开始有100只鹤，建立描述其数量变化规律的模型，并作数值计算。模型建立•记第k年沙丘鹤的数量为xk,年均增长率为r，则第k+1年鹤的数量为•xk+1=(1+r)xkk=0,1,2······•已知x0=100,在较好，中等和较差的自然环境下r=0.0194,-0.0324,和-0.0382我们利用Matlab编程，递推20年后观察沙丘鹤的数量变化情况Matlab实现•首先建立一个关于变量n，r的函数•functionx=sqh(n,r)•a=1+r;•x=100;•fork=1:n•x(k+1)=a*x(k);•end•在command窗口里调用sqh函数k=(0:20)';y1=sqh(20,0.0194);y2=sqh(20,-0.0324);y3=sqh(20,-0.0382);round([k,y1',y2',y3'])利用plot绘图观察数量变化趋势•可以用不同线型和颜色绘图•rgbcmykw分别表示红绿兰兰绿洋红黄黑白色：+o*.Xsd表示不同的线型•plot(k,y1,k,y2,k,y3)在同一坐标系下画图plot(k,y2,':')plot(k,y2,'--')plot(k,y2,'r')plot(k,y2,'y')plot(k,y2,'y',k,y1,':')plot(k,y2,k,y1,':')plot(k,y2,'oy',k,y1,':')用gtext(‘r=0.0194’),gtext(‘r=-0.0324’),gtext(‘r=-0.0382’)在图上做标记。•人工孵化是挽救濒危物种的措施之一，如果每年孵化5只鹤放入保护区，观察在中等自然条件下沙丘鹤的数量如何变化Xk+1=aXk+5,a=1+r如果我们想考察每年孵化多少只比较合适，可以令Xk+1=aXk+b,a=1+r•functionx=fhsqh(n,r,b)•a=1+r;•X=100;•Fork=1:n•X(k+1)=a*x(k)+b;•end•k=(0:20);%一个行向量•y1=(20,-0.0324,5);也是一个行向量•round([k’,y1’])对k,y1四舍五入，但是不改变变量的值plot(k,y1)ky1是行向量列向量都可以也可以观察200年的发展趋势，以及在较差条件下的发展趋势，也可以考察每年孵化数量变化的影响。一阶线性常系数差分方程的解、平衡点及其稳定性•自然环境下,b=0•人工孵化条件下•令xk=xk+1=x得差分方程的平衡点•k→∞时，xk→x,称平衡点是稳定的0kkxax10(1)kkkxaxbaa011kkaaxba1bxa1kkxaxb高阶线性常系数差分方程••如果第k+1时段变量Xk+1不仅取决于第k时段变量Xk，而且与以前时段变量有关，就要用高阶差分方程来描述一年生植物的繁殖•一年生植物春季发芽，夏天开花，秋季产种，没有腐烂，风干，被人为掠取的那些种子可以活过冬天，其中一部分能在第2年春季发芽，然后开花，产种，其中的另一部分虽未能发芽，但如又能活过一个冬天，则其中一部分可在第三年春季发芽，然后开花，产种，如此继续，一年生植物只能活1年，而近似的认为，种子最多可以活过两个冬天，试建立数学模型研究这种植物数量变化的规律，及它能一直繁殖下去的条件。模型及其求解•记一棵植物春季产种的平均数为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-2Xk=a1bcXk-1+a2b(1-a1)bcXk-2Xk=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计算了Xk=a1bcXk-1+a2b(1-a1)bcXk-2•Functionx=zwfz(x0,n,b)•C=10;a1=0.5;a2=0.25;•p=a1*b*c;q=a2*b*(1-a1)*b*c;•X1=x0;•X2=p*(x1);•fork=3:n•X(k)=p*(xk-1)+q*(xk-2);•end•K=(0:20)’;•Y1=zwfz(100,21,0.18);•Y2=zwfz(100,21,0.19);•Y3=zwfz(100,21,0,20);•Round([k,y1’,y2’,y3’])•Plot(k,y1,k,y2,’:’,k,y3,’o’),•Gtext(‘b=0.18’),gtext(‘b=0.19’),gtext(‘b=0.20’)结果分析：Xk=pXk-1+qXk-2(1)x1+px0=0(2)•对高阶差分方程可以寻求形如的解。代入(1)式得称为差分方程的特征方程。差分方程的特征根：方程(1)的解可以表为C1,c2由初始条件x0,x1确定。kkx20pq21,242ppq1122kkkxcc1,21,0()kxk•本例中，用待定系数的方法可以求出b=0.18时,c1=95.64,c2=4.36，这样实际上，植物能一直繁殖下去的条件是b0.19112(,)(0.9430,0.0430)95.64(0.9430)4.36(0.0430)kkkx1,21,()kxk1,25102b线性常系数差分方程组•汽车租赁公司的运营•一家汽车租赁公司在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.60.3•ABC•ABC•ABC假设在每个租赁期开始能把汽车都租出去，并都在租赁期末归还0.10.70.20.10.60.30.1模型及其求解•记第k个租赁期末公司在ABC市的汽车数量分别为x1(k),x2(k),x3(k)（也是第k+1个租赁期开始各个城市租出去的汽车数量），很容易写出第k+1个租赁期末公司在ABC市的汽车数量为(k=0,1,2,3···)112321233123(1)0.6()0.2()0.1()(1)0.3()0.7()0.3()(1)0.1()0.1()0.6()xkxkxkxkxkxkxkxkxkxkxkxk•用矩阵表示用matlab编程，计算x(k),观察n年以后的3个城市的汽车数量变化情况112233(1)0.60.20.1()(1)0.30.70.3()(1)0.10.10.6()xkxkxkxkxkxk•functionx=czqc(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]';•fork=1:n•x(:,k+1)=A*x(:,k);•end如果直接看10年或者20年发展趋势，可以直接在命令窗口（commondwindow）作，而不是必须编一个函数112233(1)0.60.20.1()(1)0.30.70.3()(1)0.10.10.6()xkxkxkxkxkxk•A=[0.6,0.2,0.1;0.3,0.7,0.3;0.1,0.1,0.6];•n=10;•fork=1:n•x(:,1)=[200,200,200]';•x(:,k+1)=A*x(:,k);•end•round(x)作图观察数量变化趋势012345678910120140160180200220240260280300x1(k)x2(k)x3(k)•k=0:10;•plot(k,x)，grid•gtext('x1(k)'),•gtext('x2(k)'),•gtext('x3(k)')•可以看到时间充分长以后3个城市汽车数量趋于180，300，120•可以考察这个结果与初始条件是否有关•若最开始600辆汽车都在A市，可以看到变化时间充分长以后，各城市汽车数量趋于稳定，与初始值无关直接输入x（:,1）的值即可x(:,1)=[600,0,0];round(x');plot(k,x),grid0123456789100100200300400500600按年龄分组的种群增长•野生或饲养的动物因繁殖而增加，因自然死亡和人为屠杀而减少，不同年龄动物的繁殖率，死亡率有较大差别，因此在研究某一种群数量的变化时，需要考虑年龄分组的种群增长。•将种群按年龄等间隔的分成若干个年龄组，时间也离散化为时段，给定各年龄组种群的繁殖率和死亡率，建立按年龄分组的种群增长模型，预测未来各年龄组的种群数量，并讨论时间充分长以后的变化趋势。模型及其求解•设种群按年龄等间隔的分成n个年龄组，记i=1,2,···，n,时段记作k=0,1,2···,且年龄组区间与时段长度相等(若5岁为一个年龄组，则5年为一个时段)。以雌性个体为研究对象•记在时段k第i年龄组的数量为xi(k);第i年龄组的繁殖率为bi，表示每个个体在一个时段内繁殖的数量；第i年龄组死亡率为di，表示一个时段内死亡数与总数的比，si=1-di是存活率。•注意：第k时段的第i年龄组活过来的，是第k+1时段的第i+1年龄组•Xi+1(k+1)=sixi(k)i=1，2,···,n-1,k=0,1,····•各年龄组在第k时段繁殖的数量和是第k+1时段的第1年龄组•X1(k+1)=k=0,1,····•记在时段k种群各年龄组的数量为•X(k)=[x1(k),x2(k),····,xn(k)]’1()niiibxk121121000000000nnnbbbbsLss•这样，有x(k+1)=Lx(k),k=0,1,····•给定在0时段，各年龄组的初始数量x(0)•就可以预测任意时段k,各年龄组的数量•设一种群分成5个年龄组，•繁殖率b1=0,b2=0.2,,b3=1.8,b4=0.8,b5=0.2•存活率s1=0.5,s2=0.8,s3=0.8,s4=0.1•各年龄组现有数量都是100只，•用matlab计算x(k)•b=[0,0.2,1.8,0.8,0.2];•s=diag([0.5,0.8,0.8,0.1]);•L=[b