• 検索結果がありません。

固有値問題 (1) 固有値問題に対する冪乗法 - 明治大学

N/A
N/A
Info
ダウンロード
Protected

Academic year: 2025

シェア "固有値問題 (1) 固有値問題に対する冪乗法 - 明治大学"

Copied!
7
0
0

読み込み中.... (全文を見る)

今ダウンロードする ( 7 ページ )

全文

(1)

応用数理実験

固有値問題 (1) 固有値問題に対する冪乗法

桂田 祐史 2005 年 6 月 28 日

今回は固有値問題を取り扱うための冪乗法とその周辺の種々の手法を実験してもらう。

実験に用いる行列は小さなもので良い。自分で選ぶのが面倒なら、次の行列を使ってよい1。

A =



1 1 1 1 2 2 1 2 3



.

•

選んだ行列のサイズが小さい場合は、固有多項式を計算して固有値に関する情報を調べ よ2。サイズが小さくない場合は、

Mathematica

か

(後で例をあげる)、 MATLAB (Octave)

の関数

eig()

を用いて固有値を調べよ。

•

冪乗法により絶対値最大の固有値とそれに属する固有ベクトルを求めよ。収束の速さを 調べよ。固有値を反復ベクトルの成分の比として計算した場合と、

Rayleigh

商で計算 した場合とで、収束の速さに違いが出るか？

•

逆反復法により

(連立 1

次方程式を解くのには

LU

分解を使え)、絶対値最小の固有値と それに属する固有ベクトルを求めよ。

•

シフト法の実験をせよ。すなわち絶対値最大、最小の固有値以外の固有値を見い出せ。

(例えば、上に掲げた行列 A

に対しては

0.6

に近い固有値を探してみよ。)

• (余裕がある人向け)

絶対値最大の固有値に属する固有ベクトルを求めてから、その成分

を含まない初期ベクトルを用いて冪乗法の反復を実行することにより、絶対値が二番目 に大きな固有値の固有ベクトルを求めるアルゴリズムがある。適当な文献で調べるか、

自分で考えて実験してみよ。

• (余裕がある人向け)

絶対値最大の固有値が二つ以上ある場合、冪乗法で何が起こるか観

察せよ。

1こ の 行 列 自 体 が 一 松 信「 数 値 解 析 」朝 倉 書 店 (1982) か ら 採った も の で あ る 。余 談 だ が 、こ の 本 は 基 本 的 な 数 値 計 算 ア ル ゴ リ ズ ム に 関 し て 概 観 す る に は 便 利 な 本 で あ る 。固 有 値 は λ1 = 5.048917339522305313522214407023370089866· · ·, λ2 = 0.6431041321077905561056004899786994619201· · ·, λ3= 0.3079785283699041303721851029979308815479· · ·.

2必要があればMathematica等を利用せよ。その場合のヒント: 行列式を計算するDet[],指定された次数の 単位行列を作るIdentityMatrix[],固有値を計算するEigenvalues[], 方程式を解くSolve[],方程式の解の 近似値を求めるFindRoot[],関数のグラフを描く Plot[]等の手続きがある。

(2)

Mathmatica

の使用例

¶ ³

oyabun% math

Mathematica 4.0 for Solaris

Copyright 1988-1999 Wolfram Research, Inc.

-- Motif graphics initialized -- In[1]:= A={{1,1,1},{1,2,2},{1,2,3}}

Out[1]= {{1, 1, 1}, {1, 2, 2}, {1, 2, 3}}

In[2]:= Eigenvalues[A]

Out[2]= (省略 --- 興味があればやってみて下さい) In[3]:= N[%,40]

-46 Out[3]= {5.048917339522305313522214407023369723596 + 0. 10 I,

-42

> 0.3079785283699041303721851029979308598027 + 0. 10 I, -42

> 0.6431041321077905561056004899786994166009 + 0. 10 I}

In[4]:= na=N[A,50]

Out[4]= {{1.0000000000000000000000000000000000000000000000000,

> 1.0000000000000000000000000000000000000000000000000,

> 1.0000000000000000000000000000000000000000000000000},

> {1.0000000000000000000000000000000000000000000000000,

> 2.0000000000000000000000000000000000000000000000000,

> 2.0000000000000000000000000000000000000000000000000},

> {1.0000000000000000000000000000000000000000000000000,

> 2.0000000000000000000000000000000000000000000000000,

> 3.0000000000000000000000000000000000000000000000000}}

In[5]:= Eigenvalues[na]

Out[5]= {5.0489173395223053135222144070233697235963877860565,

> 0.64310413210779055610560048997869941660087281326538,

> 0.30797852836990413037218510299793085980273940067811}

In[6]:= Eigenvectors[na]

Out[6]= {{-0.32798527760568176779603202500454990640870532112133,

> -0.59100904850610352545794571799937980844713881955014,

> -0.73697622909957824233808630700517009796156650157119},

> {0.73697622909957824233808630700517009796156650157119,

(3)

A Octave のサンプル・プログラム

グラフ化の説明もかねて

Octave

のサンプル・プログラムを以下に示す。なお、MATLAB でも同じプログラムが走ることを確認済み。

¶冪乗法 ³

format long %

表示桁数を多く取る

a=[1 1 1;1 2 2;1 2 3]

eigen=eig(a) % eig()

で固有値を求める

(カンニング)

x=ones(3,1); %

初期ベクトルを

(1,1,1)

とする

maxiter=10;

r=zeros(maxiter,1); %

近似固有値を記録するベクトルの用意

for k=1:maxiter

y=a*x;

x=y/norm(y);

r(k)=x’*(a*x); % Rayleigh

商で近似固有値を計算し、r() に記録

end

error=r-max(eigen) %

近似固有値の誤差を求める

plot(1:maxiter,abs(error)) %

反復回数と誤差をプロット

loglog(1:maxiter,abs(error)) %

両側対数目盛りでプロット

semilogy(1:maxiter,abs(error)) %

片側対数目盛りでプロット

µ ´

B MATLAB メモ

• 6701

号室の

Windows XP

マシンでは、M-file は

Z:Y

か、あるいはその下、例えば

Z:Ymatlab

に作ることを勧める。[File] メニューの

[Set Path]

サブメニューで

[Add Folder]

しておく。

• 2005

年度の情報処理教室のパソコンには、MATLAB や

Octave

がインストールされて いる。前者はライセンスの問題でつねに利用できるとは限らないが…

(4)

C 自分でやってみた

MATLAB

を使って実験してみた。非常に気楽にプログラムを書いて出来てしまう。そろそ

ろ課題を新しくして、この古い課題のプログラムをサンプルとして提示する形にするのが良い かもしれない。

C.1 冪乗法

¶ ³

% eigen1.m --- power method format long

a=[1 1 1;1 2 2;1 2 3]

eigen=eig(a) x=ones(3,1);

maxiter=10;

r=zeros(maxiter,1);

for k=1:maxiter y=a*x;

x=y/norm(y);

r(k)=x’*(a*x);

end r

error=r-max(eigen) clf

semilogy(1:maxiter,abs(error))

µ ´

1 2 3 4 5 6 7 8 9 10

10-16 10-14 10-12 10-10 10-8 10-6 10-4 10-2

(5)

C.2 Rayleigh 商の効果

¶ ³

% naive_vs_rayleigh.m format long

a=[1 1 1;1 2 2;1 2 3]

eigen=eig(a) x=ones(3,1);

maxiter=10;

r=zeros(maxiter,1);

n=zeros(maxiter,1);

for k=1:maxiter y=a*x;

x=y/norm(y);

ax=a*x;

r(k)=x’*ax;

n(k)=ax(1)/x(1);

end r n

error_r=r-max(eigen) error_n=n-max(eigen)

%subplot(1,2,1)

%semilogy(1:maxiter,abs(error_r))

%subplot(1,2,2)

%semilogy(1:maxiter,abs(error_n)) clf

semilogy(1:maxiter,abs(error_r),1:maxiter,abs(error_n))

µ ´

1 2 3 4 5 6 7 8 9 10

10-16 10-14 10-12 10-10 10-8 10-6 10-4 10-2 100

(6)

C.3 逆反復法

¶ ³

% eigen2.m --- inverse iteration method format long

a=[1 1 1;1 2 2;1 2 3]

eigen=eig(a) x=ones(3,1);

maxiter=30;

r=zeros(maxiter,1);

[L u p]=lu(a);

shift=1;

for k=1:maxiter y=u\(L\(p*x));

x=y/norm(y);

r(k)=x’*(a*x);

end r

error=r-min(eigen) clf

semilogy(1:maxiter,abs(error))

µ ´

0 5 10 15 20 25 30

10-16 10-14 10-12 10-10 10-8 10-6 10-4 10-2 100

(7)

C.4 シフト法

¶ ³

% eigen3.m --- shift method format long

a=[1 1 1;1 2 2;1 2 3]

eigen=eig(a) x=ones(3,1);

maxiter=10;

r=zeros(maxiter,1);

shift=0.6;

% we cannot use the factorization of a, can we?

[L u p]=lu(a-shift*eye(3,3));

for k=1:maxiter y=u\(L\(p*x));

x=y/norm(y);

r(k)=x’*(a*x);

end r

eigen=sort(eigen);

error=r-eigen(2) clf

semilogy(1:maxiter,abs(error))

µ ´

1 2 3 4 5 6 7 8 9 10

10-16 10-14 10-12 10-10 10-8 10-6 10-4 10-2

参照

今ダウンロードする ( PDF - 7 ページ - 91.52 KB )

関連したドキュメント

有限オートマトンでの計算可能性問題 • 言語 A ⊂ Σ∗ に対し

後置記法の演算式のスタックを用いた計算 逆ポーランド電卓 • 数値 =⇒ push • 演算子 =⇒ 被演算子を所定の個数だけ pop −→ 演算を施し、結果を push • 入力終了 =⇒ pop −→ スタックが丁度空になったらその値が答え 問 : 後置記法逆ポーランド記法の式に対し スタックを用いて値を計算する

ハバードモデルの超大規模固有値問題に対する地球シミュレータでの並列計算法 (数値解析と新しい情報技術)

MPI (Message Passing Interface) を利用した並列化を行なう. 4.1 方法 1

大規模一般化固有値問題の解法とその並列化 (微分方程式の数値解法と線形計算)

Nagoya University 1 はじめに 本論文では一般化固有値問題 $Ax=\lambda Bx$ において , $\sim$ 行列 $A,$

ホモトピー法に基づく2固有値問題の数値解法の収束性解析(科学技術における数値計算の理論と応用)

数値解法 前節において、

の未定乗数法 説明と練習・基本問題

という関係があるとする。変数 を固定したとき、この関係式は に対する方程式にな る。もしこの方程式が唯一の解 を持つとき、実数

並列分散計算によるBMI固有値問題解法

Ninf は，クライアントからサーバに対して計算の 実行を依頼する RPC 型の Grid 計算システムであ る．Ninf による Grid 計算では，クライアント計算

複素関数・同演習第 14 - 明治大学

Mathematica の勧め 現象数理学科 Mac にインストールされている Mathematica は、数式処理 系と呼ばれるソフトウェアである。プログラミング言語処理系の一種でもある が、多くのプログラミング言語 例えば C, Python, MATLAB, … は、数値計算 はできても数式の計算はできない。

行列の固有値問題

この符号の変化数は特別な注意をせずに 数値的に安定して計算できる。つまり絶対値が 非常に小さくて、符号の判別がつきにくい場合も、「符号の変化数」そのものは疑いがなく計 算できる。例えば pk−1x pkx pk+1x + 絶対値小 − または pk−1x pkx pk+1x − 絶対値小 + において pkx が正であっても負であっても 0