2019年差分方程及微分方程数值解.doc

东南大学《数学实验》报告
实验内容：差分方程及微分方程数值解 一 实验目的
熟悉迭代法及微分方程数值方法 二 预备知识
（1）了解差分方程稳定性、周期分解、混沌等相关知识 （2）了解欧拉方法、龙格-库特方法。

三 实验内容与要求
（一）V olterra 方程数值解
方程
0,,,,,>⎪⎪⎩⎪⎪⎨
⎧+-=-=d c b a dxy cx dt
dy bxy ax dt
dx
其中
a=1,b=0.1,c=0.5,d=0.04
.2)0(,25)0(==y x
命令与结果
在函数编辑器中输入：
function dxdt = euler( t,x ) dxdt=[ x(1)*(1-0.1*x(2)) x(2)*(-0.5+0.02*x(1))]; end
四阶龙格-库塔公式：
在命令窗口中输入： tspan=[0 15]; x0=[25;2];
[t,x]=ode45(@euler,tspan,x0); plot(t,x(:,1),'r-','LineWidth',0.5); hold on;
plot(t,x(:,2),'g-','LineWidth',0.5); hold on;
（1） 相图
四阶龙格-库塔公式：
axis([0 15 0 125])
legend('x(1)','x(2)')
grid on
title('The Numerical Solution Of
Calculable way of fourth rank
Rounge-kutt')
plot(x(:,1),x(:,2))
欧拉公式：
欧拉方法：
在命令窗口中输入：
tspan=[0 15];
x0=[25;2];
[t,x]=ode23(@euler,tspan,x0);
plot(t,x(:,1),'r-','LineWidth',0.5);
hold on;
plot(t,x(:,2),'g-','LineWidth',0.5);
hold on;
axis([0 15 0 125])
legend('x(1)','x(2)')
grid on
title('The Numerical Solution Of
Euler Equation')
>> plot(x(:,1),x(:,2))
（2）轨线图
四/五阶龙格-库塔公式：
欧拉公式：
（二）差分阻滞增长模型
在t 时刻单位时间内的人口数量的变化量仅仅与此时的人口数量x 有关(等于右边的值),
)
1()(N
x
rx t x -=
其中的r 表示人口的固有增长率,N 表示能容纳的最大人口数.
经过 单位时间，即一个繁殖周期的种群数量的增长量仅仅与前一个时期的种群数量有关。

于是模型为
即
这是一个一阶非线性差分方程.对于给定的初值,我们可以从这个递推公式运用计算机很容易
地计算出一些xk.这是在计算机出现以后的一个新的特点.但是我们更关心的是当时间趋于
无穷时,即k 趋于无穷时, xk 的极限如何，即差分方程平衡点的稳定性问题。

求得差分方程(2)的平衡点为0和(b-1)/b.
它们分别对应于差分方程(1)的平衡点0和N. 我们将这个差分方程(2)在平衡点附近展开,有 k
k bx x x ==+1:0处
注意到b=1+r>1,平衡点0是不稳定的.
*)
)(2(*)(*:1
1*21x x b x x b x x b x k k k --+--=--=+
略去高阶项得
*)
)(2(*1x x b x x k k --=-+
因此当|2-b|<1时,平衡点x*＝ (b-1)/b 是稳定的. |2-b|>1即 b>3时,平衡点x*是不稳定的.
（1）稳定情况（b=1.6,2.2）
b=1.6 b=2.2
(1) ,2,1,0),1(1 =-=-+k N y ry y y k
k k k )(1 '=+-+=+ ,2,1,0),)1(1()1(1k N r ry y r y k
k k 为两边乘以则上式记令))1((,1.)1(N r r
r b N r ry x k
k ++=+=
(2)
)1(1k k k x bx x -=+
b=1.6,在1到3范围之间 稳定
命令与结果（~()n x n 图）
function x=block(b,x0) b=1.6;
if nargin<2 x0=0.5; end x=x0;
n=1:10000; for i=2:10000
x(i)=b*x(i-1)*(1-x(i-1)); end
plot(n,x,'.'); end
结果~()n x n 图：
b=2.2,在1到3范围之间 稳定
命令与结果（~()n x n 图）
function x=block(b,x0) b=2.2;
if nargin<2 x0=0.5; end x=x0;
n=1:10000; for i=2:10000
x(i)=b*x(i-1)*(1-x(i-1)); end
plot(n,x,'.'); end
结果~()n x n 图：
（2）不稳定情况
b=3.3 b=3.63
b=3.3在1到3范围之间
不稳定
命令与结果（~()
n x n图）b=3.63在1到3范围之间
不稳定
命令与结果（~()
n x n图）
（3）倍周期分解
b=3.3b=3.5
（1）求稳定平衡点
可得两个平衡点：0.4794 0.8236 程序：
function x=block(b,x0)
b=3.3;
if nargin<2
x0=0.5;
end
x=x0;
n=1:100;
for i=2:100
x(i)=b*x(i-1)*(1-x(i-1));
end
（2）数值结果
ans =
Columns 1 through 16
0.5000 0.8250 0.4764 0.8232 0.4804 0.8237 0.4792 0.8236 0.4795 0.8236 0.4794 0.8236
（1）求稳定平衡点
四个平衡点：0.5009 0.8750 0.3828 0.8269
程序：
function x=block(b,x0)
b=3.5;
if nargin<2
x0=0.5;
end
x=x0;
n=1:100;
for i=2:100
x(i)=b*x(i-1)*(1-x(i-1));
end
（2）数值结果
ans =
Columns 1 through 16
0.5000 0.8750 0.3828 0.8269 0.5009 0.8750 0.3828 0.8269
0.4794 0.8236 0.4794 0.8236 Columns 17 through 32
0.4794 0.8236 0.4794 0.8236 0.4794 0.8236 0.4794 0.8236 0.4794 0.8236 0.4794 0.8236 0.4794 0.8236 0.4794 0.8236
……0.5009 0.8750 0.3828 0.8269 0.5009 0.8750 0.3828 0.8269
Columns 17 through 32
0.5009 0.8750 0.3828 0.8269 0.5009 0.8750 0.3828 0.8269 0.5009 0.8750 0.3828 0.8269 0.5009 0.8750 0.3828 0.8269
……。