差分与微分方程

差分方程
理论：
1.一阶差分方程
k k k x x x -=∆+1….刻画该变量
形如)(k k x f x =∆或)(1k k x F x =+称为一阶差分方程；
2.二阶差分方程
k k k k k k x x x x x x +-=∆-∆=∆+++12122
形如),(12k k k x x f x +=∆称为二阶差分方程
3.平衡点和稳定性
如果*lim x x k k =+∞
→即平衡点 渐近稳定：存在*x 的某个邻域U ，对任意的U x ∈0,虽然*0x x ≠，但*lim x x k k =+∞
→ 4.应用及软件实现：
一阶线性常系数差分方程，)1.1(,......2,1,0,)1(1=+=+k x r x k k 其中r 为常数，有3种方式计算k 时段的增长率
前差公式：k
k
k x x x -+1 中点公式：k
k k x x x 21
1-+- 后差公式：k k k x x x 1
--
其中中点公式的精度最高
)1.1(的解为等比数列,......2,1,0,)1(0=+=k r x x k k
若0≠r ，则仅有平衡点0=x 。

稳定当且仅当1|1|<+r
下面选取参数r 和初始值0x ，按)1.1(迭代，绘图观察其解的长期行为 详见程序
r=[0.09;0.09;-0.1;-0.1;-1.9;-1.9;-2.09;-2.09];
x=[15;-15;85;-85;85;-85;15;-15];
一阶线性常系数非齐次差分方程)2.1(,......2,1,0,)1(1=++=+k b x r x k k
若0=r 则为等差数列0,0,1,2,......k x x kb k =+=；
若0≠r ，则r
b r r b x x k k -++=)1)((0 引入 r
b x y k k +=则.0,1,2.....k )1()1(01=+=+=+k k k r y y r y 可得此时平衡点r
b x -=稳定当且仅当02-<<r 实例：Florida 沙丘鹤属于濒危物种，生态学家估计它在较好的自然环境下，年平均增长率仅为 1.94%，而在中等及较差自然环境下年平均增长率仅为-3.24%和-3.82%，即它逐渐减少，假设在某自然保护区内开始时有100只沙丘鹤，请建立数学模型，描述其数量变化规律，并作数值计算。

人工孵化是挽救濒危物种的措施之一，如果每年人工孵化5只沙丘鹤放入该保护区，问在3中自然条件下沙丘鹤的数量将会如何变化？
r=[0.09;0.09;-0.1;-0.1;-1.9;-1.9;-2.09;-2.09];
x=[15;-15;85;-85;85;-85;15;-15];
for n=1:20
x(:,n+1)=(1+r).*(x(:,n));
end
s{1}='x0>0,r>0';
s{2}='1';
s{3}='1';
s{4}='1';
s{5}='1';
s{6}='1';
s{7}='1';
s{8}='1';
for k=1:8
subplot(4,2,k),plot(0:20,x(k,:),'k+')
axis([-1,21,-100,100])
xlabel(s{k})
end
r=[0.0194,-0.0324,-0.0382];
x=[100,100,100];
for k=1:20
x(k+1,:)=x(k,:).*(1+r);
end
disp('yan bian')
disp(' year good median bad')
disp([(0:20)',round(x)])
plot(0:20,x(:,1),'k^',0:20,x(:,2),'ko',0:20,x(:,3),'kv')
思考问题：
贷款有两种还款方式：等额本金和等额本息
设有房款60万，首付12万，贷款48万，期限20年，月利率0.5%，问两种还款方式最终差多少钱？
二、某种山猫在较好中等以及较差的自然环境下，年平均增长率分别为：1.68%，0.055%和-4.5%，假设开始时有100只山猫，按以下情况分别讨论山猫数量逐年变化的过程及趋势
1.三种自然环境下25年的变化情况，结果要列表并图示；
2.如果每年捕获3只，山猫数量将会如何变化？会灭绝吗？如果每年捕获一只呢？
3.在较差自然环境下，如果要使山猫数量稳定在60只左右，每年要人工繁殖多少只？
二阶线性常系数齐次差分方程：
形式：......2,1,012=+=++k bx ax x k k k 其中04,0,,2≠+≠∈b a b R b a
其特征方程为02=--b a λλ
所以有2个互异的根(一对实根或共轭复根) 21,λλ，因为0,0,021≠≠∴≠λλb 于是一般解为1122,0,1,2......k k k x c c k λλ=+=
其中21,c c 为任意常数，给定10,x x 得到
⎩⎨⎧+=+=2
2111210λλc c x c c x 得到唯一解 1≠+b a 则仅有平衡点0=x
可以证明当同时1||,1||21<<λλ时平衡点0=x 是稳定的
斐波那契数列：在一年之初把一对一雌一雄新生的兔子放入围栏，从第二个月开始，母兔每月生出一对一雌一雄的小兔，每对新生的兔子也从它们第二个月开始，每月生出一对一雌一雄的小兔，求一年后围栏内有多少只兔子？第二年后呢？
从理论和软件两方面分析
酵母培养物的增长
问题提出：
下表数据是从测量酵母培养物增长的实验收集而来，请建立数学模型，模拟酵母培养物的增长过程
K=[0:18];
xk=[9.6,18.3,29.0,47.2,71.1,119.1,174.6,257.3,350.7,441.0,513.3,559.7,594.8,629.4,640.8,651.1,655.9,661.8];
问题分析：
首先绘制xk 关于k 的散点图
极限k k x +∞
→lim 可能存在 再计算k x ∆并绘图
前图是k x ∆关于k,后图是关于xk ，看出后者呈现明显的二次曲线关系 利用前差公式计算k k
k k x x x r -=+1 (1)
前图是rk 关于k ，后图是rk 关于xk 的散点图以及直线拟合
易见rk 关于xk 近似线性递减关系
0,0)/1(>>-=N r N x r r k k (2)
将（1）代入（2），即建立了离散阻滞增长模型
)1(1N
x rx x x k k k k -+=+ 符号说明：
K ：小时；
Xk ：酵母培养物在k 小时的生物量；
rk ：生物量在k 小时的增长率；
r ：生物的固有增长率；
N ：生物的最大容纳量
建立模型：
1(1)k k k k x x x rx N
+=+- 模型求解和模型检验：
需要参数r,N 以及初始值x0，
下面介绍两种方法：
1.由rk 关于xk 的散点图以及直线拟合，设rk=a(2)+a(1)xk
当rk=0时，认为xk 达到最大值)
1()2(a a - 当xk=0，认为是固有增长率)2(a r =
-40-20
20
40
效果不是特别理想，误差有些大
2.利用非线性拟合
sse2 =
1.3535e+003 0200
400
600
024681012141618
-40-20
20
40
练习：一、给出美国人口的预报(2010-2020)单位：百万
Year=[1790:10:2000];
Renkou=[3.9,5.3,7.2,9.6,12.9,17.1,23.2,31.4,38.6,50.2,62.9,76,92,106.5,123.2,131.7,150.7,179.3,204.0,226.5,251.4,281.4];
两种方式预报：
1.线性递减函数0,0)/1(>>-=N r N x r r k k
.2.指数衰减函数31
2a e a r k x a k +=- 哪种方式更精确？
酵母培养程序：
01一些绘图命令
k=[0:18];
xk=[9.6,18.3,29.0,47.2,71.1,119.1,174.6,257.3,350.7,441.0,513.3,559.7,594.8,629.4,640.8,651.1,655.9,659.6,661.8];
plot(k,xk,'k+')
for i=2:19
delxk(i)=xk(i)-xk(i-1);
rk(i)=(xk(i)-xk(i-1))./xk(i-1);
end
figure(2)
subplot(1,2,1),plot(k,delxk,'k+')
subplot(1,2,2),plot(xk,delxk,'k+')
hold on
a=polyfit(xk,delxk,2);
subplot(1,2,2),plot(polyval(a,0:750))
axis([0,750,0,100]);
figure(3)
subplot(1,2,1),plot(k,rk,'k+')
subplot(1,2,2),plot(xk,rk,'k+')
hold on
a=polyfit(xk,rk,1);
subplot(1,2,2),plot(polyval(a,0:750))
axis([0,750,0,1]);
02经验初值
t=[0:18];
x=[9.6,18.3,29.0,47.2,71.1,119.1,174.6,257.3,350.7,441.0,513.3,559.7,594.8,629.4,640.8,651.1,655.9,659.6,661.8];
r=(x(2:19)-x(1:18))./x(1:18);
a1=polyfit(x(1:18),r,1);
r1=a1(2)
N1=-a1(2)./a1(1)
x1=x(1);
for k=1:18
x1(k+1)=x1(k)+r1*x1(k)*(1-x1(k)/N1);
end
resdl=x-x1;
ssel=sum(resdl.^2)
axis([-1,19,0,670]);
subplot(2,1,1),plot(t,x,'k+',t,x1,'ks')
legend('观测值','模拟值')
subplot(2,1,2),plot(t,resdl,'k+',[-1,19],[0,0],'k:')
axis([-1,19,-40,40]);
03 非线性拟合
clear;
clc;
t=[0:18];
x=[9.6,18.3,29.0,47.2,71.1,119.1,174.6,257.3,350.7,441.0,513.3,559.7,594.8,629.4,64 0.8,651.1,655.9,659.6,661.8];
r=(x(2:19)-x(1:18))./x(1:18);
[a2,se2,rj]=nlinfit(t,x,@fun_jiaomu,[0.5,660,9.6])
sse2=sum(se2.^2)
subplot(2,1,1),plot(t,x,'k+',t,fun_jiaomu(a2,t),'ks')
axis([-1,19,0,670]);
legend('观测值','模拟值')
subplot(2,1,2),plot(t,se2,'k+',[-1,19],[0,0],'k')
axis([-1,19,-40,40]);
函数
%b(1)=r,b(2)=N,b(3)=x_0
function y=fun_jiaomu(b,x)
y=zeros(size(x));
y(1)=b(3);
for k=2:length(x)
y(k)=y(k-1)+b(1).*y(k-1).*(1-y(k-1)./b(2));
end。