matlab迭代。混沌,原因,分叉 Microsoft Word 文档
- 1、下载文档前请自行甄别文档内容的完整性,平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益,请联系客服反馈,我们会尽快为您处理(人工客服工作时间:9:00-18:30)。
问题与实验3: 一元线性迭代的收敛性条件怎样表述? 关于迭代法收敛性的两个判别条件: a 、充分必要条件是:矩阵M 的谱半径
(){}1,..,2,1max
<==n i M i
i
λ
ρ()
b 、充分条件是:矩阵M 的某个算子范数M
<1。
问题与实验4: 在本例中,12
<M
,这时迭代序列是收敛的,就本例
或选择别的例子,按
12
<M
和12
≥M
构造不同的迭代法,通过实验
和比较,并给出你对实验结果的解释(如关于收敛性、收敛速度等),当然这需要你首先知道矩阵范数的概念,并且对它有比较好的理解。
设x 是方程组(5)的解,{}m
x 是迭代法(6)生成的任一序列,
因为
f Mx x +=,f Mx
x m
m +=+1
()()()
022
1x x M
x x M
x x M x x m
m m m -
==-
=-=--- ,
设D = diag (a 11, a 22, …, a nn ),将AX = b 改写为: AX = (D – (D - A )) x = b DX = (D - A ) x + b
X = (I – D -1A ) x + D -1b
记 B = I – D -1A F = D -1 b 则迭代格式的向量表示为
F BX X
k k +=+)
()
1( B称为雅克比迭代矩阵。
由此可知要判断X 是否收敛只需看M 的谱半径是否小于1,既有一
其中I 是单位i 矩阵,D 是提取A 的对角线上的元素。
下判断条件:充要条件:(1) (){}1,..,2,1max
<==n i M i
i
λ
ρ.(2)
充分
条件是:矩阵M 的某个算子范数
M
<1.并且我们知道当M 越小的时
候其收敛的速度越快。
并且还可以知道当初始值越接近精确解时收敛速度越快。
这是由于迭代的公式所定的。
下面来看另一个例子:X1+2X2+-2X3=1 X1+X2+X3=1 2X1+2X2+X3=1
雅可比法的迭代矩阵:A=[1 2 -2; 1 1 1; 2 2 1;] b=[1;1;1;] D=diag(diag(A)); LU=D-A; M=D\LU;
p=max(abs(eig(M))) f=D\b; x=[]; z=[];
x(:,1)=eye(3,1); N=200000000; for i=1:N;
if norm(A*x(:,i)-b)<1e-010;
m=i;break
else
x(:,i+1)=M*x(:,i)+f;
z=x(:,i+1)
end
end
m
e=norm(A*z-b)
plot([1:length(x)],x)
title('JACOBI ITERATION OF LINEAR EQUATIONS')
A =1 2 -2
1 1 1
2 2 1
b = 1 1 1
p = 5.8106e-006(谱范数)可以看出是收敛的
z=1 0 -1
z =-1 1 -1
z =-33 1
m =4(迭代的次数)
e =0(误差的估计)
图像是:
1 1.5
2 2.5
3 3.54
-3
-2
-1
1
2
3
JACOBI ITERATION OF LINEAR EQUATIONS
A=[9 -1 -1; -1 8 0; -1 0 9;] b=[7;7;8;] D=diag(diag(A)); LU=D-A; M=D\LU;
p=max(abs(eig(M))) f=D\b; x=[]; z=[];
x(:,1)=eye(3,1);
N=200000000;
for i=1:N;
if norm(A*x(:,i)-b)<1e-010;
m=i;break
else
x(:,i+1)=M*x(:,i)+f;
z=x(:,i+1)
end
end
m
e=norm(A*z-b)
plot([1:length(x)],x)
title('JACOBI ITERATION OF LINEAR EQUATIONS')
A = 9 -1 -1
-1 8 0
-1 0 9
b =7 7 8
p =0.1620(谱范数)可以看出是收敛,m =16迭代的次数。
(雅可比法的迭代)e = 1.7316e-011其误差范围。
z = 1.0000 1.0000 1.0000 其图像是:p = 5.8106e-006(谱范数) m =4(迭代的次数)
2
4
6
8
10
12
14
16
00.10.20.30.40.50.60.70.80.9
1JACOBI ITERATION OF LINEAR EQUATIONS
2.高斯――赛得尔(Gauss-Seidel )迭代法
显然,如果迭代收敛,
)
1(+k i
x 应该比
)
(k i
x 更接近于原方程的解*
i
x (i =
1, 2,…, n ),因此在迭代过程中及时地以
)
1(+k i
x 代替
)
(k i
x (i = 1, 2,…,
n-1),可望收到更好的效果。
这样(4.2)式可写成:
()()()⎪⎪⎪
⎪⎩
⎪⎪⎪
⎪⎨
⎧---=---=---=+--++++++)
1(11,)1(22)1(11)1()
(2)(323)1(121222)1(2)
(1)(313)(212111
)1(1111k n n n k n k n n nn
k n k n
n k k k k n n k k k x a x a x a b a x x a x a x a b a x x a x a x a b a x
(4.5)
(4.5)式可简写成
⎪⎪⎭
⎫
⎝
⎛
-
-=
∑∑+=+-=+)(1
)1(11
)
1(1k j
n
i j ij
k j
i j ij
i ii k i
x
a
x
a
b a x (i = 1, 2,…, n )
此为G -S 迭代格式。
G -S 迭代格式的矩阵表示:
Jacob 方法的迭代格式为
F
Bx
x
k k +=+)
()
1(
其中 B = I-D -1A (I 是单位矩阵) F = D -1b 将B 写成
B = L + U L 是下三角阵(表示是从
X1从到XN 依次求解),U 是上三角阵 从而 F
Ux Lx
x
k k k ++=++)
()
1()
1(
F
Ux
x L I k k +=-+)
()
1()(
F
L I Ux
L I x
k k 1
)
(1
)
1()()(--+-+-=
(4.6)
此时迭代矩阵为 (I - L )-1U
高斯――赛得尔(Gauss-Seidel )迭代法编程代码:
function x=Gauss_Seidel(A,b,x0,tol) if (nargin==2) x0=ones(size(b)); tol=1e-6;
elseif (nargin==3) tol=1e-6; else
sprintf('USAGE:Gauss_Seidel(A,b,x0,tol)') end
D=diag(diag(A)); U=triu(A,1);
L=tril(A,-1);
G=-(D+L)\U;
d1=(D+L)\b;
x=G*x0+d1;
n=1;
while norm(x-x0)>=tol
x0=x;
x=G*x0+d1;
n=n+1;
end
n
运行结果是:A=[9 -1 -1;
-1 8 0;
-1 0 9;]
b=[7;7;8;]
x=gauss_seidel(A,b)
A = 9 -1 -1
-1 8 0
-1 0 9
b =7 7 8
n = 1(运行的次数)收敛的数度是比较快的。
x = 1 1 1
A= [4.3999 0.6686 -1.6041 0.5287 -1.0106; 0.6900 5.1908 0.2573 0.2193 0.6145; 0.8156 -1.2025 4.0565 -0.9219 0.5077; 0.7120 -0.0198 1.4151 3.1707 1.6924; 1.2902 -0.1567 -0.8051 -0.0592 5.5913] b=[ -0.4326; -1.6656;0.1253;0.2877;-1.1465] x=gauss_seidel(A,b)
A =4.3999 0.6686 -1.6041 0.5287 -1.0106 0.6900 5.1908 0.2573 0.2193 0.6145 0.8156 -1.2025 4.0565 -0.9219 0.5077 0.7120 -0.0198 1.4151 3.1707 1.6924 1.2902 -0.1567 -0.8051 -0.0592 5.5913 b = -0.4321 -1.6656 0.1253 0.2877 -1.1465
n = 14(运行的次数。
)可以看出用高斯――赛得尔(Gauss-Seidel )
的收敛速度比雅可比迭代法(42次)快。
051015202530354045
-0.5
0.5
1
JACOBI IT ERAT ION OF LINEAR EQUAT IONS
x1
x2x3
x4x5
x = -0.1070 -0.2950 0.0322 0.1957 -0.1819
问题与实验5:对迭代(
)2
1
a
x a x
n n -
-=+,通过进一步的实验和观
察寻找更多的k —循环的情况,并借助于必要的理论分析,讨论次迭代出现k —循环的原因。
1.迭代法 (不动点迭代法)
迭代法是求解一元非线性方程 0)(=x f 的主要方法。
其做法是(类似线性代数方程组的迭代法):将方程0)(=x f 改写为等价方程 )(x x ϕ=
这时,方程)(x x ϑ=成为“隐式”形式,除非ϑ是x 的线性函数,否则不能直接算出它的根。
对此,从某个初始值0
x 开始,对
应)(x x ϕ=构造迭代公式(格式)
)
,2,.1,0)((1 ==+k x x k k ϕ
这就可确定序列{})
,2,1,0( =k
x
k
(因x 是单变量,所以迭代标记
k
写为下标即可)。
如果{}k
x 有极限
*
lim x
x k k =∞
→
并且由 )
(1
k k x x
ϕ=+式两边取极限,可知*
x 就是方程)(x x ϕ=的
根,这时称迭代公式)
(1
k k x x ϕ=+是收敛的。
实际求解中,只有
对收敛的)
(1
k k x x
ϕ=+才有意义,但不能做无穷多步,因此,按
精度要求,取某个迭代值k
x 作为方程)(x x ϕ=的根的近似值,这种求根的方法称为迭代法,式中 )(x ϕ称为迭代函数。
这里,由于)
(*
*
x x
ϕ=,直观看来,即解*
x 经ϕ计算后)(*
x ϕ仍
等于*
x (岂不是没有动),因此,称*
x 是函数ϕ的一个不动
⋅⋅⋅= ,,,210n
点,公式)
,2,.1,0)((1
==+k x x
k k ϕ也称为不动点迭代公式,非线性
方程的迭代法也称为不动点迭代法。
二、求不动点的迭代格式 定义6 设.若存在点,使得有
g(
)=
.则称
为映射
的一个不动点.(迭代一次
时), 不动点是曲线
与直线
的交
点.(表示不断的迭代。
)
.如果已知在上连续,
且能证明存在极限,则由
,可知
的极限即为在
中的一个不动点.然而上述
不一定收敛,下图为
时四
种迭代过程的示意图,其中(a)、(b)收敛,(c),(D )则发散.
(a)
(b)
其中(a)、( b)收敛, (c).(d)则发散.
(c) 发散. (d) 发散.
从几何直观上看,(c )
发散的原因可能是
在不动点
附近的变化率太大所致.下面给出的压缩映射原理将证实这一直观猜想.
三、压缩映射原理 定义7 设
,且满足:
(i) ;
(ii)
使
这时称为上的一个压缩映射.
1.一些基本概念
设)(x f 是一个定义在实数域上的实值函数,如果存在*
x ,使得*
*
)(x x f =,则称*
x 为)(x f 的不动点。
如果所有附近的点在迭代过程中都趋向于某个不动点,则称该不动点为吸引点,或称为稳定点;如果所有附近的点在迭代过程中都远离他而去,则称该不动点为排斥点,或称为不稳定点。
如果
1
13221)(,)(,,)(,)(a a f a a f a a f a a f k k k ====- ,
且k j a a j ,,3,2,1 =≠,则k a a a ,,,21 构成一个k 循环。
1
a 称为k 周期点,k
a a a ,,,21 称为一个k 周期轨道。
0 <λ ≤ 4 , 0 < x ≤ 1。
在上面典型的非线性迭代方程中,还发现有“倍周期分叉”现象:▲ 当 1< λ < 3时, 迭代的归宿是一个确定的数ξ 。
例如:λ=2.4时, x n+1=x n =7/12(n →∞),周期为1。
▲ 当 λ≥3时, 迭代出现多个确定的数值ξ 。
(λ = 3 时,曲线开始分叉)
例如:λ = 3.2 时,一个λ值对应的有两个ξ 值,
)(n
n
n x x x -=+1 1
λ⋅⋅⋅= ,,,210n
即其归宿轮流取两个值:x n +2 = x n , 0.7955 → 0.5130, 周期为2。
λ = 3.5 时,与一个λ 值对应的有4个ξ 值,
即其归宿轮流取四个值x n +4=x n 0.3828→0.8269
↑ ↓ 周期为4。
0.8750←0.5009 ▲ 当 3.570≤λ≤4 时,最后归宿可取无穷多值, 即出现混沌现象,周期为∞。
通过倍周期分叉走向混沌的道路,这是目前已知的一种典型方式,如下图所示:
引例 关于迭代
(
)x a x a
n n +=--
12
(7)
实验与讨论。
容易求得迭代(7)有两个不动点:0*
=x 和12*
-=a x ,图4给出
了a
=2
及迭代初值分别取01
00
.=x
和900
.=x
时,迭代收敛于1
2*
-=a x 的情形.(是一循环) 使用程序(ITERA1:) a=2; x=[];y=[];
x(1)=input('please input the firse value of iteration x(1)=:');%x(1)=0.01
y(1)=input('please input the firse value of iteration y(1)=:');%y(1)=0.9
n=input('please input the total number of iteration n=:');%n=150 for i=1:n-1;
x(i+1)=a-(x(i)-sqrt(a))^2; y(i+1)=a-(y(i)-sqrt(a))^2; end
subplot(2,1,1) plot(1:n,x,'bo')
title('ITERATION FOR SOLVING EQUATION') subplot(2,1,2) plot(1:n,y,'r*')
title('ITERATION FOR SOLVING EQUATION')
50
100
150
00.511.52ITERATION FOR SOLVING EQUATION
050100150
0.5
1
1.5
2ITERATION FOR SOLVING EQUATION
图5给出了a
3时,出现4—循环的情形:
使用程序(ITERA2:) a=3; x=[];
x(1)=input('please input the firse value of iteration x(1)=:');%2.5
n=input('please input the total number of iteration n=:');%n=200 for i=1:n-1;
x(i+1)=a-(x(i)-sqrt(a))^2; end
plot(1:n,x,'r*')
title('ITERATION FOR SOLVING EQUATION')
020406080100120140160180200
1.4
1.61.82
2.22.42.62.83ITERATION FOR SOLVING EQUATION
图 5
function fencha
%分叉图,x0为起点,a 为参数
x(1)=2.5; for a=3:0.01:4 for i=1:200
x(i+1)=a-(x(i)-sqrt(a))^2; if i>30
plot(a,x(i))
hold on end end end
3
3.1
3.2
3.3
3.4
3.5
3.6
3.7
3.8
3.9
4
00.5
1
1.5
2
2.5
3
3.5
4
问题与实验6:对迭代(
)2
1
a
x a x
n n -
-=+,通过进一步的实验和观察确
定可以出现混沌的参数a 的取值范围,混沌现象的出现是否与迭代初值有关?
1.一些基本概念
设)(x f 是一个定义在实数域上的实值函数,如果存在*
x ,使得*
*
)(x x f =,则称*
x 为)(x f 的不动点。
如果所有附近的点在迭代过程中都趋向于某个不动点,则称该不动点为吸引点,或称为稳定点;如果所有附近的点在迭代过程中都远离他而去,则称该不动点为排斥点,或称为不稳定点。
如果
1
13221)(,)(,,)(,)(a a f a a f a a f a a f k k k ====- ,
且k j a a j
,,3,2,1
=≠,则k
a a a ,,,2
1
构成一个k 循环。
1
a 称为k 周期点,k
a a a ,,,21 称为一个k 周期轨道。
0 <λ ≤ 4 , 0 < x ≤ 1。
在上面典型的非线性迭代方程中,还发现有“倍周期分叉”现象:▲ 当 1< λ < 3时, 迭代的归宿是一个确定的数ξ 。
例如:λ=2.4时, x n+1=x n =7/12(n →∞),周期为1。
▲ 当 λ≥3时, 迭代出现多个确定的数值ξ 。
(λ = 3 时,曲线开始分叉)
例如:λ = 3.2 时,一个λ值对应的有两个ξ 值, 即其归宿轮流取两个值:x n +2 = x n , 0.7955 → 0.5130, 周期为2。
λ = 3.5 时,与一个λ 值对应的有4个ξ 值,
即其归宿轮流取四个值x n +4=x n 0.3828→0.8269
↑ ↓ 周期为4。
0.8750←0.5009 ▲ 当 3.570≤λ≤4 时,最后归宿可取无穷多值, 即出现混沌现象,周期为∞。
通过倍周期分叉走向混沌的道路,这是目前已知的一种典型方式,如下图所示:
)(n
n
n x x x -=+1 1
λ⋅⋅⋅= ,,,210n
3
3.1
3.2
3.3
3.4
3.5
3.6
3.7
3.8
3.9
4
00.5
1
1.5
2
2.5
3
3.5
4
可以看出K 循环及混沌的出现是与a 有关的,与初始值无关。
x(n+1)=a -(x(n)-√2)^2的图像其中a=3的分叉图
a=3.5; x=[];
x(1)=input('please input the firse value of iteration x(1)=:');%1.5
n=input('please input the total number of iteration n=:');%n=500 for i=1:n-1;
x(i+1)=a-(x(i)-sqrt(a))^2; end
plot(1:n,x,'r*')
title('ITERATION FOR SOLVING EQUATION')
please input the firse value of iteration x(1)=:0.9(初始的值) please input the total number of iteration n=:200
050100150200250300
0.5
1
1.5
2
2.5
3
3.5
IT ERAT ION FOR SOLVING EQUAT ION
050100150200250300
0.5
1
1.5
2
2.5
3
3.5
ITERATION FOR SOLVING EQUATION
a=3.5x(1)=0.9(初
始的值)n=:200
a=3.5x(1)=3.4(初始
的值)n=300在收敛的时候,其收敛与否与初值无关。
:
50
100
150
01234ITERATION FOR SOLVING EQUATION
050100150
1234ITERATION FOR SOLVING EQUATION
50
100
150
200
250
300
01234ITERATION FOR SOLVING EQUATION
50
100
150
200
250
300
01234ITERATION FOR SOLVING EQUATION
x=[];y=[];
x(1)=input('please input the firse value of iteration x(1)=:');%x(1)=0.50001 y(1)=input('please input the firse value of iteration y(1)=:');%y(1)=0.9 n=input('please input the total number of iteration n=:');%n=150
下面看混沌的情况:图7给出了a
=4
时,出现混沌现象的情形,迭代初值分别取
900.=x 50001
00.=x 。
混沌现象的产生与初值的选取无关。
3.7出现K 循环x(1)=3 y(1)=0.9 n=300 a=3.7
for i=1:n-1;
x(i+1)=a-(x(i)-sqrt(a))^2;
y(i+1)=a-(y(i)-sqrt(a))^2;
end
subplot(2,1,1)
plot(1:n,x,'bo')
title('ITERATION FOR SOLVING EQUATION') subplot(2,1,2)
plot(1:n,y,'r*')
title('ITERATION FOR SOLVING EQUATION')。