正面教材分享 http://blog.sciencenet.cn/u/wdlang 70%的以色列人是无神论者,不过他们都相信上帝给了他们那块土地。这个世界经不起思考

博文

计算方法8:幂法

已有 6045 次阅读 2016-12-27 20:32 |个人分类:计算方法|系统分类:教学心得

公告:各位同学,期末考试在1月10号。有不懂的问题尽早来问。

幂法用来计算一个矩阵的模最大的本征值及其本征向量。

其操作极其简单。一个n乘n矩阵,其实是一个从某个线性空间到其自身的线性变换在某个基下的表示。因此,一个矩阵A最天然的功能就是对一个矢量v作用,将之变成另外一个矢量w。对这个新的w,我们又可以用A作用。如此迭代,最终矢量会收敛到与模最大的那个本征值相对应的本征向量。

这里的情况有点类似卓别林在电影摩登时代里的做法。他有个扳手。一个扳手最自然的作用就是拧一个螺丝。所以当他手里有个扳手的时候,他就到处找螺丝,最后都找到在等公共汽车的路人的裤子铝扣上去。

我们手里有个矩阵,其最自然的作用是乘以一个矢量。所以我们能做的就是找个矢量,对其作用作用再作用。

幂法简单直观,但是其潜在的推广的空间很大。幂法的不足是它只保留当前的矢量,而弃之前的矢量于不顾,好比猴子掰玉米。要是我们保留之前的矢量,我们就能从中提取有关矩阵的更多的信息。这就是krylov的思想。从矢量v开始产生的一系列矢量v, Av, A^2v, A^3v,...张成一个子空间,此即所谓krylov空间。幂法很单纯地取最后产生的那个矢量作为要求的本征矢量的近似。但是一个更聪明的想法是取所有这些矢量的一个线性叠加作为要求的本征矢量的近似。这便是严格对角化的lanczos方法的做法。

以下程序可以演示在幂法下矢量如何收敛到要求的矢量。

% power method

clear all; close all; clc;


delta = 7;

H = [1+delta, 1; 1, 1+delta];

v0 = [1;0];

N = 25;


vlist = zeros(2, N);

vlist(:,1) = v0;


h1= plot(nan);

hold on

plot(v0(1), v0(2),'o')

xlim([0,1])

ylim([0,1])

axis square

mov(1) = getframe;


for s = 2: N

   temp = H* vlist(:,s-1);

   temp = temp/norm(temp);

   vlist(:,s) = temp;

   

   set(h1,'XData',vlist(1,1:s),'YData',vlist(2,1:s))

   set(h1,'Marker','o')

   drawnow

   pause(1)

   mov(s) = getframe;

end


str=strcat('delta=',num2str(delta),'.avi');

movie2avi(mov, str, 'compression', 'None','fps',2);




https://wap.sciencenet.cn/blog-100379-1023709.html

上一篇:计算方法7:利用幅角原理确定一个整函数在给定区域内的零点数
下一篇:PRA审稿意见
收藏 IP: 103.214.68.*| 热度|

2 姬扬 徐令予

该博文允许注册用户评论 请点击登录 评论 (2 个评论)

数据加载中...

Archiver|手机版|科学网 ( 京ICP备07017567号-12 )

GMT+8, 2024-4-25 09:55

Powered by ScienceNet.cn

Copyright © 2007- 中国科学报社

返回顶部