Ｃ言語による数値計算プログラミング演習

Size: px

Start display at page:

Download "Ｃ言語による数値計算プログラミング演習"

さあしゃかんざとばる
7 years ago
Views:

1 5. 行列の固有値問題 n n 正方行列 A に対する n 個の固有値 λ i (i=1,,,n) と対応する固有ベクトル u i は次式を満たす Au = λ u i i i a11 a1 L a1 n u1i a1 a a n u i A =, ui = M O M M an 1 an L ann uni これらはまとめて, つぎのように書ける 5.1 ヤコビ法 = Λ, = [ u1 u u n ], Λ = AU U U L 実数対称行列 A の固有値と対応する固有ベクトルをすべて求める方法であり, 古典的ではあるが, 行列サイズが 10 以下ならば現在でも使用されている直交行列 P を用いて A の相似変換 ( すなわち直交変換 )Ã=P T AP を行い, 対角行列に収束させていくと, 対角要素が固有値となる ( 相似変換に関して固有値は不変 ) λ1 0 λ O 0 λ n アルゴリズム A:= 与えられた行列 U:=I( 単位行列 ) for m=1,, { 収束条件を満たす ( すべての非対角要素がεより小さくなる ) まで繰り返す for i=1 to n-1 {A の対角要素を除いた上三角部分を走査 for j=i+1 to n if a ij ε then { 非対角要素で絶対値がεより大きい要素 a ij を零にする変換 θ:=(1/)cot -1 ((a ii -a jj )/(a ij )) A:=P T AP {A の直交変換 ( 非対角要素 a ij を零にする ) U:=UP {U の変換 n n- n-1 n

2 A の直交変換 Ã=P T AP はつぎのように書き表せ, cotθ = (a ii -a jj )/(a ij ) ã ii = (1/)(a ii +a jj ) + (1/)(a ii -a jj )cosθ + a ij sinθ ã jj = (1/)(a ii +a jj ) - (1/)(a ii -a jj )cosθ - a ij sinθ ã ij = ã ji = 0 ã ik = ã ki = a ik cosθ + a jk sinθ (k=1,,,n, ただし k i,j) ã jk = ã kj = -a ik sinθ + a jk cosθ (k=1,,,n, ただし k i,j) ã kl = a kl (k,l=1,,,n, ただし k,l i,j) U の変換 Ũ=UP はつぎのようになる ũ ki = u ki cosθ + u kj sinθ (k=1,,,n) ũ kj = -u ki sinθ + u kj cosθ (k=1,,,n) ũ kl = u kl (k,l=1,,,n, ただし l i,j) 計算誤差を少なくするための注意 : 実際の数値計算においては,θ を求めることなく, 代数式から cosθ, sinθ, cosθ, sinθ を求めるほうが, 計算時間, 計算誤差とも少なくできる一般に,cotα( あるいは tanα=cot -1 α) から cosα, sinα を求める場合, cotα 1 のとき sin α= 1, cosα = cotα sinα 1+ cot α cotα 1 のとき tanα = 1 1, cos α = 1, sinα = tanα cosα cotα 1+ tan α により計算すると,1 よりはるかに大きい数の二乗を行わずにすむのでオーバーフローを避けることができ, かつ 0 による割算を避けることができる cotθ から cosθ, sinθ を求めるときはこの関係を用いるさらに cosθ, sinθ を求めるときには半角公式を用いるが, 桁落ちを避けるため倍角公式も併用して, cos θ 0 1 のとき cos θ = (1 + cos θ), sinθ = sin θ (cos θ ) cosθ 0 1 のとき sin θ = (1 cos θ), cosθ = sin θ (sin θ ) により計算する実用に耐える数値計算では, このような細部にまで注意が払われていることに留意しようプログラム /* solve eigen-value problem for a symmetric matrix */ /* by Jacobi's method */ #include <stdio.h> #include <math.h> #define NMAX 0 int jacobi(float); float a[nmax][nmax], u[nmax][nmax];

3 int n; void main(void) { int mcum,i,j; float eps; printf("order of a Symmetric matrix: n? "); scanf("%d",&n); printf("input Lower triangle of matrix n"); printf("? A(%d, j), j=1,%d : ", i+1, i+1); for(j=0;j<=i;j++){ scanf("%g",&a[i][j]); a[j][i]=a[i][j]; printf(" nsymmetric matrix: A n"); for(j=0;j<n;j++) printf(" %15.6e", a[i][j]); printf(" n"); for(j=0;j<i;j++) u[i][j]=u[j][i]=0; u[i][i]=1; mcum= 0; /* initial set for cumulative number of transformations */ while(printf(" ntolerance : eps? "), scanf("%g",&eps)!=eof){ mcum += jacobi(eps); printf("cumulative no. of transformations=%d n", mcum); printf(" n"); for(i=0;i<n;i++) printf(" %15.6e", a[i][i]); printf(" n n"); for(j=0;j<n;j++) printf(" %15.6e", u[i][j]); printf(" n"); int jacobi(float eps) { int i,j,k,m,mcycle; float p,q,t,s,c,c,s,r; m=0; /* initial set for transformation counter (total) */ 3

4 while(1){ mcycle=0; /* initial set for transformation counter (one cycle) */ for(i=0;i<n-1;i++) for(j=i+1;j<n;j++) if( fabs(a[i][j])>=eps ){ m++; mcycle++; /* increment of transformation counters */ p= (a[i][i]-a[j][j])/; q=a[i][j]; if( fabs(p)<fabs(q) ){ t=p/q; /* t=cot( theta) */ s= 1/sqrt(1+t*t); if(q<0) s= -s; c= t*s; else{ t=q/p; /* t=tan( theta) */ c= 1/sqrt(1+t*t); if(p<0) c= -c; s= t*c; if(c>0){ c= sqrt( (1+c)/ ); s= s/c/; else{ s= sqrt( (1-c)/ ); c= s/s/; r= (a[i][i]+a[j][j])/; a[i][i]= r + p*c + q*s; a[j][j]= r - p*c - q*s; a[i][j]= a[j][i]=0; for(k=0;k<n;k++) if( (k!=i)&&(k!=j) ){ p= a[i][k]; q=a[j][k]; a[i][k]= a[k][i]= p*c+q*s; a[j][k]= a[k][j]= -p*s+q*c; for(k=0;k<n;k++){ p= u[k][i]; q= u[k][j]; u[k][i]= p*c+q*s; u[k][j]=-p*s+q*c; if(mcycle==0) return(m); /* return no. of transformations (total) */ 実行例教科書章末問題の質点ばね系の固有値問題をヤコビ法により解く (1) 運動方程式はつぎのように行列ベクトル表示される 4

5 M&& s = Ks s1 m1 0 0 k1+ k k 0 s = s, M = 0 m 0, K = k k + k k 3 3 s m3 0 k3 k3 () 系は固有振動数 ω の調和運動をしているとすれば, 解は振幅を s = xe iωt と表され, これを運動方程式に代入すると, 一般固有値問題 K x = ω M x となるので, 結局, 標準固有値問題 Ax 1 = λx ( ここに, A= M K, λ = ω ) を得る m 1 =m =m 3 =m, k 1 =k =k 3 =k とすると, 行列 A は 1 0 k A = 1 1 m と表されるいま k/m=1 として, 標準固有値問題を解く実行開始 Order of a Symmetric matrix: n? 3 Input Lower triangle of matrix? A( 1, j), j=1, 1 :? A(, j), j=1, : -1? A( 3, j), j=1, 3 : Symmetric matrix: A e e e e e e e e e+00 Tolerance : eps? 0.1 cumulative no. of transformations= e e e e e e e e e e e e-01 Tolerance : eps? 0.01 cumulative no. of transformations= e e e e e e e e e e e e-01 Tolerance : eps? x T = [ x1, x, x] として

6 cumulative no. of transformations= e e e e e e e e e e e e-01 Tolerance : eps? 1.0e-6 cumulative no. of transformations= e e e e e e e e e e e e-01 Tolerance : eps? ^Z (control と Z を同時に押す ) (Enter Key を押す ) 入力終了コードの入力おしまい参考文献葉子 : 数値計算の基礎解法と誤差, コロナ社 (007) 森口繁一, 伊理正夫, 武市正人編 :C による算法痛論, 東京大学出版会 (000) Heath, Michael T.: Scientific Computing, An Introductory Survey, McGraw-Hill(00) 6

Ｃ言語による数値計算プログラミング演習

Ｃ言語による数値計算プログラミング演習 8. 数値積分任意の区間 [,b] における f() の定積分 b () I = f ( ) d の値は, つぎのように n 点の関数値の和により近似的に与えられる () In = Ak f ( k) n k = このとき, k を分点,A k を重みという 8. ニュートンコーツ ( 複合型 ) 積分公式積分区間 [,b] を等分割して n 個の分点をとり, 被積分関数 f() を n- 次ラグランジュ補間多項式で近似して得られる積分公式を