a = 1 b = 2 printf ("%12.10f\n" , bisection(a, b)) def bisection(a, b) { while (true) { // 区間 (a, b) の中点 c = (a + b) / 2 c = (a + b) / 2 printf("%12.10f\t%13.10f\n", c, c - Math.sqrt(2)) fc = f(c) if (Math.abs(fc) < 0.0000000001) break if (fc < 0) { // f(c) < 0 であれば, 解は区間 (c, b) の中に存在 a = c } else { // f(c) > 0 であれば, 解は区間 (a, c) の中に存在 b = c } } c } def f(x) { x * x - 2 }
a = 1 b = 2 printf ("%12.10f\n" , falseposition(a, b)) def falseposition(a, b) { while (true) { // 点 (a,f(a)) と 点 (b,f(b)) を結ぶ直線と x軸の交点 c = (a * f(b) - b * f(a)) / (f(b) - f(a)) printf("%12.10f\t%13.10f\n", c, c - Math.sqrt(2)) fc = f(c) if (Math.abs(fc) < 0.0000000001) break if (fc < 0) { // f(c) < 0 であれば, 解は区間 (c, b) の中に存在 a = c } else { // f(c) > 0 であれば, 解は区間 (a, c) の中に存在 b = c } } c } def f(x) { x * x - 2 }
x = 1 printf ("%12.10f\n" , iterative(x)) def iterative(x0) { while (true) { x1 = g(x0) printf("%12.10f\t%13.10f\n", x1, x1 - Math.sqrt(2)) if (Math.abs(x1 - x0) < 0.0000000001) break x0 = x1 } x1 } def g(x) { (x / 2) + (1 / x) }
x = 2 printf ("%12.10f\n" , newton(x)) def newton(x0) { while (true) { x1 = x0 - (f0(x0) / f1(x0)) printf("%12.10f\t%13.10f\n", x1, x1 - Math.sqrt(2)) if (Math.abs(x1 - x0) < 0.0000000001) break x0 = x1 } x1 } def f0(x) { x * x - 2 } def f1(x) { 2 * x }
x = 2 printf ("%12.10f\n" , bailey(x)) def bailey(x0) { while (true) { x1 = x0 - (f0(x0) / (f1(x0) - (f0(x0) * f2(x0) / (2 * f1(x0))))) printf("%12.10f\t%13.10f\n", x1, x1 - Math.sqrt(2)) if (Math.abs(x1 - x0) < 0.0000000001) break x0 = x1 } x1 } def f0(x) { x * x - 2 } def f1(x) { 2 * x } def f2(x) { 2 }
x0 = 1 x1 = 2 printf ("%12.10f\n" , secant(x0, x1)) def secant(x0, x1) { while (true) { x2 = x1 - f(x1) * (x1 - x0) / (f(x1) - f(x0)) printf("%12.10f\t%13.10f\n", x2, x2 - Math.sqrt(2)) if (Math.abs(x2 - x1) < 0.0000000001) break x0 = x1 x1 = x2 } x2 } def f(x) { x * x - 2 }
N = 3 def a = [[ 9.0, 2.0, 1.0, 1.0], [2.0, 8.0, -2.0, 1.0], [-1.0, -2.0, 7.0, -2.0], [1.0, -1.0, -2.0, 6.0]] as double[][] def b = [20.0, 16.0, 8.0, 17.0] as double[] def c = [ 0.0, 0.0, 0.0, 0.0] as double[] // ヤコビの反復法 jacobi(a,b,c) println("X") disp_vector(c) // ヤコビの反復法 def jacobi(a, b, x0) { def x1 = [ 0.0, 0.0, 0.0, 0.0] as double[] while (true) { def finish = true for (i in 0..N ) { x1[i] = 0 for (j in 0..N ) { if (j != i) x1[i] += a[i][j] * x0[j] } x1[i] = (b[i] - x1[i]) / a[i][i] if (Math.abs(x1[i] - x0[i]) > 0.0000000001) finish = false } for (i in 0..N ) x0[i] = x1[i] if (finish) return disp_vector(x0) } } // 1次元配列を表示 def disp_vector(row) { for (col in row) printf ("%14.10f\t" , col) println() }
N = 3 def a = [[ 9.0, 2.0, 1.0, 1.0], [2.0, 8.0, -2.0, 1.0], [-1.0, -2.0, 7.0, -2.0], [1.0, -1.0, -2.0, 6.0]] as double[][] def b = [20.0, 16.0, 8.0, 17.0] as double[] def c = [ 0.0, 0.0, 0.0, 0.0] as double[] // ガウス・ザイデル法 gauss(a,b,c) println("X") disp_vector(c) // ガウス・ザイデル法 def gauss(a, b, x0) { def x1 = [ 0.0, 0.0, 0.0, 0.0] as double[] while (true) { def finish = true for (i in 0..N ) { x1[i] = 0 for (j in 0..N ) { if (j != i) x1[i] += a[i][j] * x0[j] } x1[i] = (b[i] - x1[i]) / a[i][i] if (Math.abs(x1[i] - x0[i]) > 0.0000000001) finish = false x0[i] = x1[i] } if (finish) return disp_vector(x0) } } // 1次元配列を表示 def disp_vector(row) { for (col in row) { printf ("%14.10f\t" , col) } println() }
N = 3 def a = [[ 9.0, 2.0, 1.0, 1.0], [2.0, 8.0, -2.0, 1.0], [-1.0, -2.0, 7.0, -2.0], [1.0, -1.0, -2.0, 6.0]] as double[][] def b = [20.0, 16.0, 8.0, 17.0] as double[] def c = [ 0.0, 0.0, 0.0, 0.0] as double[] // ピボット選択 pivoting(a,b) println("pivoting") println("A") disp_matrix(a) println("B") disp_vector(b) println() // 前進消去 forward_elimination(a,b) println("forward elimination") println("A") disp_matrix(a) println("B") disp_vector(b) println() // 後退代入 backward_substitution(a,b) println("X") disp_vector(b) // 前進消去 def forward_elimination(a, b) { for (pivot in 0..(N - 1)) { for (row in (pivot + 1)..N) { s = a[row][pivot] / a[pivot][pivot] for (col in pivot..N) a[row][col] -= a[pivot][col] * s b[row] -= b[pivot] * s } } } // 後退代入 def backward_substitution(a, b) { for (row = N; row >= 0; row--) { for (col = N; col > row; col--) b[row] -= a[row][col] * b[col] b[row] /= a[row][row] } } // ピボット選択 def pivoting(a, b) { for (pivot in 0..N) { // 各列で 一番値が大きい行を 探す max_row = pivot max_val = 0 for (row in pivot..N) { if (Math.abs(a[row][pivot]) > max_val) { // 一番値が大きい行 max_val = Math.abs(a[row][pivot]) max_row = row } } // 一番値が大きい行と入れ替え if (max_row != pivot) { tmp = 0 for (col in 0..N) { tmp = a[max_row][col] a[max_row][col] = a[pivot][col] a[pivot][col] = tmp } tmp = b[max_row] b[max_row] = b[pivot] b[pivot] = tmp } } } // 1次元配列を表示 def disp_matrix(matrix) { for (row in matrix) { for (col in row) { printf ("%14.10f\t" , col) } println() } } // 2次元配列を表示 def disp_vector(row) { for (col in row) { printf ("%14.10f\t" , col) } println() }
N = 3 def a = [[-1.0, -2.0, 7.0, -2.0], [1.0, -1.0, -2.0, 6.0], [ 9.0, 2.0, 1.0, 1.0], [2.0, 8.0, -2.0, 1.0]] as double[][] def b = [ 8.0, 17.0, 20.0, 16.0] as double[] def c = [ 0.0, 0.0, 0.0, 0.0] as double[] // ピボット選択 pivoting(a,b) println("pivoting") println("A") disp_matrix(a) println("B") disp_vector(b) println() // 前進消去 forward_elimination(a,b) println("forward elimination") println("A") disp_matrix(a) println("B") disp_vector(b) println() // 後退代入 backward_substitution(a,b) println("X") disp_vector(b) // 前進消去 def forward_elimination(a, b) { for (pivot in 0..N) { for (row in 0..N) { if (row == pivot) continue s = a[row][pivot] / a[pivot][pivot] for (col in pivot..N) a[row][col] -= a[pivot][col] * s b[row] -= b[pivot] * s } } } // 後退代入 def backward_substitution(a, b) { for (pivot in 0..N) b[pivot] /= a[pivot][pivot] } // ピボット選択 def pivoting(a, b) { for (pivot in 0..N) { // 各列で 一番値が大きい行を 探す max_row = pivot max_val = 0 for (row in pivot..N) { if (Math.abs(a[row][pivot]) > max_val) { // 一番値が大きい行 max_val = Math.abs(a[row][pivot]) max_row = row } } // 一番値が大きい行と入れ替え if (max_row != pivot) { tmp = 0 for (col in 0..N) { tmp = a[max_row][col] a[max_row][col] = a[pivot][col] a[pivot][col] = tmp } tmp = b[max_row] b[max_row] = b[pivot] b[pivot] = tmp } } } // 1次元配列を表示 def disp_matrix(matrix) { for (row in matrix) { for (col in row) { printf ("%14.10f\t" , col) } println() } } // 2次元配列を表示 def disp_vector(row) { for (col in row) { printf ("%14.10f\t" , col) } println() }
N = 3 def a = [[-1.0, -2.0, 7.0, -2.0], [1.0, -1.0, -2.0, 6.0], [ 9.0, 2.0, 1.0, 1.0], [2.0, 8.0, -2.0, 1.0]] as double[][] def b = [ 8.0, 17.0, 20.0, 16.0] as double[] def c = [ 0.0, 0.0, 0.0, 0.0] as double[] // ピボット選択 pivoting(a,b) println("pivoting") println("A") disp_matrix(a) println("B") disp_vector(b) println() // 前進消去 forward_elimination(a,b) println("LU") disp_matrix(a) // Ly=b から y を求める (前進代入) y = [0.0,0.0,0.0,0.0] forward_substitution(a,b,y) println "Y" disp_vector(y) // Ux=y から x を求める (後退代入) x = [0.0,0.0,0.0,0.0] backward_substitution(a,y,x) println("X") disp_vector(x) // ピボット選択 def pivoting(a, b) { for (pivot in 0..N) { // 各列で 一番値が大きい行を 探す max_row = pivot max_val = 0 for (row in pivot..N) { if (Math.abs(a[row][pivot]) > max_val) { // 一番値が大きい行 max_val = Math.abs(a[row][pivot]) max_row = row } } // 一番値が大きい行と入れ替え if (max_row != pivot) { tmp = 0 for (col in 0..N) { tmp = a[max_row][col] a[max_row][col] = a[pivot][col] a[pivot][col] = tmp } tmp = b[max_row] b[max_row] = b[pivot] b[pivot] = tmp } } } // 前進消去 def forward_elimination(a, b) { for (pivot in 0..(N - 1)) { for (row in (pivot + 1)..N) { s = a[row][pivot] / a[pivot][pivot] for (col in pivot..N) a[row][col] -= a[pivot][col] * s // これが 上三角行列 a[row][pivot] = s // これが 下三角行列 // b[row] -= b[pivot] * s // この値は変更しない } } } // Ly=b から y を求める (前進代入) def forward_substitution(a, b, y) { for (row in 0..N) { for (col in 0..row) b[row] -= a[row][col] * y[col] y[row] = b[row] } } // Ux=y から x を求める (後退代入) def backward_substitution(a, y, x) { for (row = N; row >= 0; row--) { for (col = N; col > row; col--) y[row] -= a[row][col] * x[col] x[row] = y[row] / a[row][row] } } // 1次元配列を表示 def disp_matrix(matrix) { for (row in matrix) { for (col in row) { printf ("%14.10f\t" , col) } println() } } // 2次元配列を表示 def disp_vector(row) { for (col in row) { printf ("%14.10f\t" , col) } println() }
N = 3 def a = [[5.0,2.0,3.0,4.0],[2.0,10.0,6.0,7.0],[3.0,6.0,15.0,9.0],[4.0,7.0,9.0,20.0]] as double[][] def b = [34.0,68.0,96.0,125.0] as double[] println("A") disp_matrix(a) println("B") disp_vector(b) println() // 前進消去 forward_elimination(a,b) println("LL^T") disp_matrix(a) println() // Ly=b から y を求める (前進代入) y = [0.0,0.0,0.0,0.0] forward_substitution(a,b,y) println "Y" disp_vector(y) // L^Tx=y から x を求める (後退代入) x = [0.0,0.0,0.0,0.0] backward_substitution(a,y,x) println("X") disp_vector(x) // 前進消去 def forward_elimination(a, b) { for (pivot in 0..N) { s = 0 if (pivot > 0) for (col in 0..(pivot-1)) s += a[pivot][col] * a[pivot][col] // ここで根号の中が負の値になると計算できない! a[pivot][pivot] = Math.sqrt(a[pivot][pivot] - s) if (pivot >= N) continue for (row in (pivot+1)..N) { s = 0 if (pivot > 0) for (col in 0..(pivot-1)) s += a[row][col] * a[pivot][col] a[row][pivot] = (a[row][pivot] - s) / a[pivot][pivot] a[pivot][row] = a[row][pivot] } } } // 前進代入 def forward_substitution(a, b, y) { for (row in 0..N) { for (col in 0..row) b[row] -= a[row][col] * y[col] y[row] = b[row] / a[row][row] } } // 後退代入 def backward_substitution(a, y, x) { for (row = N; row >= 0; row--) { for (col = N; col > row; col--) y[row] -= a[row][col] * x[col] x[row] = y[row] / a[row][row] } } // 1次元配列を表示 def disp_matrix(matrix) { for (row in matrix) { for (col in row) { printf ("%14.10f\t" , col) } println() } } // 2次元配列を表示 def disp_vector(row) { for (col in row) { printf ("%14.10f\t" , col) } println() }
N = 3 def a = [[5.0,2.0,3.0,4.0],[2.0,10.0,6.0,7.0],[3.0,6.0,15.0,9.0],[4.0,7.0,9.0,20.0]] as double[][] def b = [34.0,68.0,96.0,125.0] as double[] println("A") disp_matrix(a) println("B") disp_vector(b) println() // 前進消去 forward_elimination(a,b) println("LDL^T") disp_matrix(a) println() // Ly=b から y を求める (前進代入) y = [0.0,0.0,0.0,0.0] forward_substitution(a,b,y) println "Y" disp_vector(y) // DL^Tx=y から x を求める (後退代入) x = [0.0,0.0,0.0,0.0] backward_substitution(a,y,x) println("X") disp_vector(x) // 前進消去 def forward_elimination(a, b) { for (pivot in 0..N) { // pivot < k の場合 s = 0 if (pivot > 0) for (col in 0..(pivot-1)) { s = a[pivot][col] if (col > 0) for (k in 0..(col-1)) s -= a[pivot][k] * a[col][k] * a[k][k] a[pivot][col] = s / a[col][col] a[col][pivot] = a[pivot][col] } // pivot == k の場合 s = a[pivot][pivot] if (pivot > 0) for (k in 0..(pivot-1)) s -= a[pivot][k] * a[pivot][k] * a[k][k] a[pivot][pivot] = s } } // 前進代入 def forward_substitution(a, b, y) { for (row in 0..N) { for (col in 0..row) b[row] -= a[row][col] * y[col] y[row] = b[row] } } // 後退代入 def backward_substitution(a, y, x) { for (row = N; row >= 0; row--) { for (col = N; col > row; col--) y[row] -= a[row][col] * a[row][row] * x[col] x[row] = y[row] / a[row][row] } } // 1次元配列を表示 def disp_matrix(matrix) { for (row in matrix) { for (col in row) { printf ("%14.10f\t" , col) } println() } } // 2次元配列を表示 def disp_vector(row) { for (col in row) { printf ("%14.10f\t" , col) } println() }
N = 3 // ベキ乗法で最大固有値を求める main() def main() { def a = [[5.0, 4.0, 1.0, 1.0], [4.0, 5.0, 1.0, 1.0], [1.0, 1.0, 4.0, 2.0], [1.0, 1.0, 2.0, 4.0]] as double[][] def x = [1.0 ,0.0 ,0.0 ,0.0] as double[] // ベキ乗法 lambda = power(a, x) println() println("eigenvalue") printf("%14.10f\n" , lambda) println("eigenvector") disp_vector(x) } // ベキ乗法 def power(a, x0) { lambda = 0.0 // 正規化 (ベクトル x0 の長さを1にする) normarize(x0) e0 = 0.0 for (i in 0..N) e0 += x0[i] for (k in 1..200) { // 1次元配列を表示 printf("%3d\t" , k) disp_vector(x0) // 行列の積 x1 = A × x0 x1 = [0.0 ,0.0 ,0.0 ,0.0] for (i in 0..N) for (j in 0..N) x1[i] += a[i][j] * x0[j] // 内積 p0 = 0.0 p1 = 0.0 for (i in 0..N) { p0 += x1[i] * x1[i] p1 += x1[i] * x0[i] } // 固有値 lambda = p0 / p1 // 正規化 (ベクトル x1 の長さを1にする) normarize(x1) // 収束判定 e1 = 0.0 for (i in 0..N) e1 += x1[i] if (Math.abs(e0 - e1) < 0.00000000001) break for (i in 0..N) x0[i] = x1[i] e0 = e1 } lambda } // 1次元配列を表示 def disp_vector(row) { for (col in row) { printf("%14.10f\t" , col) } println() } // 正規化 (ベクトルの長さを1にする) def normarize(x) { s = 0.0 for (i in 0..N) s += x[i] * x[i] s = Math.sqrt(s) for (i in 0..N) x[i] /= s }
N = 3 // 逆ベキ乗法で最小固有値を求める main() def main() { def a = [[5.0, 4.0, 1.0, 1.0], [4.0, 5.0, 1.0, 1.0], [1.0, 1.0, 4.0, 2.0], [1.0, 1.0, 2.0, 4.0]] as double[][] def x = [1.0 ,0.0 ,0.0 ,0.0] as double[] // LU分解 forward_elimination(a) // 逆ベキ乗法 lambda = inverse(a, x) println() println("eigenvalue") printf("%14.10f\n" , lambda) println("eigenvector") disp_vector(x) } // 逆ベキ乗法 def inverse(a, x0) { lambda = 0.0 // 正規化 (ベクトル x0 の長さを1にする) normarize(x0) e0 = 0.0 for (i in 0..N) e0 += x0[i] for (k in 1..200) { // 1次元配列を表示 printf("%3d\t" , k) disp_vector(x0) // Ly=b から y を求める (前進代入) b = [0.0,0.0,0.0,0.0] y = [0.0,0.0,0.0,0.0] for (i in 0..N) b[i] = x0[i] forward_substitution(a,y,b) // Ux=y から x を求める (後退代入) x1 = [0.0,0.0,0.0,0.0] backward_substitution(a,x1,y) // 内積 p0 = 0.0 p1 = 0.0 for (i in 0..N) { p0 += x1[i] * x1[i] p1 += x1[i] * x0[i] } // 固有値 lambda = p1 / p0 // 正規化 (ベクトル x1 の長さを1にする) normarize(x1) // 収束判定 e1 = 0.0 for (i in 0..N) e1 += x1[i] if (Math.abs(e0 - e1) < 0.00000000001) break for (i in 0..N) x0[i] = x1[i] e0 = e1 } lambda } // 前進消去 def forward_elimination(a) { for (pivot in 0..(N - 1)) { for (row in (pivot + 1)..N) { s = a[row][pivot] / a[pivot][pivot] for (col in pivot..N) a[row][col] -= a[pivot][col] * s // これが 上三角行列 a[row][pivot] = s // これが 下三角行列 } } } // Ly=b から y を求める (前進代入) def forward_substitution(a, y, b) { for (row in 0..N) { for (col in 0..row) b[row] -= a[row][col] * y[col] y[row] = b[row] } } // Ux=y から x を求める (後退代入) def backward_substitution(a, x, y) { for (row = N; row >= 0; row--) { for (col = N; col > row; col--) y[row] -= a[row][col] * x[col] x[row] = y[row] / a[row][row] } } // 1次元配列を表示 def disp_vector(row) { for (col in row) { printf("%14.10f\t" , col) } println() } // 正規化 (ベクトルの長さを1にする) def normarize(x) { s = 0.0 for (i in 0..N) s += x[i] * x[i] s = Math.sqrt(s) for (i in 0..N) x[i] /= s }
N = 3 // LR分解で固有値を求める main() def main() { def a = [[5.0, 4.0, 1.0, 1.0], [4.0, 5.0, 1.0, 1.0], [1.0, 1.0, 4.0, 2.0], [1.0, 1.0, 2.0, 4.0]] as double[][] def l = [[0.0 ,0.0 ,0.0 ,0.0], [0.0 ,0.0 ,0.0 ,0.0], [0.0 ,0.0 ,0.0 ,0.0], [0.0 ,0.0 ,0.0 ,0.0]] as double[][] def u = [[0.0 ,0.0 ,0.0 ,0.0], [0.0 ,0.0 ,0.0 ,0.0], [0.0 ,0.0 ,0.0 ,0.0], [0.0 ,0.0 ,0.0 ,0.0]] as double[][] for (k in 1..200) { // LU分解 decomp(a, l, u) // 行列の積 multiply(u, l, a) // 対角要素を表示 printf("%3d\t" , k) disp_eigenvalue(a) // 収束判定 e = 0.0 for (i in 1..N) for (j in 0..(i - 1)) e += Math.abs(a[i][j]) if (e < 0.00000000001) break } println() println("eigenvalue") disp_eigenvalue(a) } // LU分解 def decomp(a, l, u) { for (i in 0..N) { for (j in 0..N) { l[i][j] = 0.0 u[i][j] = 0.0 } } l[0][0] = 1.0 for (j in 0..N) u[0][j] = a[0][j] for (i in 1..N) { u[i][0] = 0.0 l[0][i] = 0.0 l[i][0] = a[i][0] / u[0][0] } for (i in 1..N) { l[i][i] = 1.0 t = a[i][i] for (k in 0..i) t -= l[i][k] * u[k][i] u[i][i] = t if (i < N) { for (j in (i + 1)..N) { u[j][i] = 0.0 l[i][j] = 0.0 t = a[j][i] for (k in 0..i) t -= l[j][k] * u[k][i] l[j][i] = t / u[i][i] t = a[i][j] for (k in 0..i) t -= l[i][k] * u[k][j] u[i][j] = t } } } } // 行列の積 def multiply(a, b, c) { for (i in 0..N) { for (j in 0..N) { s = 0.0 for (k in 0..N) s += a[i][k] * b[k][j] c[i][j] = s } } } // 対角要素を表示 def disp_eigenvalue(a) { for (i in 0..N) printf("%14.10f\t" , a[i][i]) println() }
N = 3 // QR分解で固有値を求める main() def main() { def a = [[5.0, 4.0, 1.0, 1.0], [4.0, 5.0, 1.0, 1.0], [1.0, 1.0, 4.0, 2.0], [1.0, 1.0, 2.0, 4.0]] as double[][] def q = [[0.0 ,0.0 ,0.0 ,0.0], [0.0 ,0.0 ,0.0 ,0.0], [0.0 ,0.0 ,0.0 ,0.0], [0.0 ,0.0 ,0.0 ,0.0]] as double[][] def r = [[0.0 ,0.0 ,0.0 ,0.0], [0.0 ,0.0 ,0.0 ,0.0], [0.0 ,0.0 ,0.0 ,0.0], [0.0 ,0.0 ,0.0 ,0.0]] as double[][] for (k in 1..200) { // QR分解 decomp(a, q, r) // 行列の積 multiply(r, q, a) // 対角要素を表示 printf("%3d\t" , k) disp_eigenvalue(a) // 収束判定 e = 0.0 for (i in 1..N) for (j in 0..(i - 1)) e += Math.abs(a[i][j]) if (e < 0.00000000001) break } println() println("eigenvalue") disp_eigenvalue(a) } // QR分解 def decomp(a, q, r) { x = [0.0 ,0.0 ,0.0 ,0.0] as double[] for (k in 0..N) { for (i in 0..N) x[i] = a[i][k] if (k > 0) { for (j in 0..(k - 1)) { t = 0.0 for (i in 0..N) t += a[i][k] * q[i][j] r[j][k] = t r[k][j] = 0.0 for (i in 0..N) x[i] -= t * q[i][j] } } s = 0.0 for (i in 0..N) s += x[i] * x[i] r[k][k] = Math.sqrt(s) for (i in 0..N) q[i][k] = x[i] / r[k][k] } } // 行列の積 def multiply(a, b, c) { for (i in 0..N) { for (j in 0..N) { s = 0.0 for (k in 0..N) s += a[i][k] * b[k][j] c[i][j] = s } } } // 対角要素を表示 def disp_eigenvalue(a) { for (i in 0..N) printf("%14.10f\t" , a[i][i]) println() }
N = 3 // ヤコビ法で固有値を求める main() def main() { def a = [[5.0, 4.0, 1.0, 1.0], [4.0, 5.0, 1.0, 1.0], [1.0, 1.0, 4.0, 2.0], [1.0, 1.0, 2.0, 4.0]] as double[][] def v = [[1.0 ,0.0 ,0.0 ,0.0], [0.0 ,1.0 ,0.0 ,0.0], [0.0 ,0.0 ,1.0 ,0.0], [0.0 ,0.0 ,0.0 ,1.0]] as double[][] // ヤコビ法 jacobi(a, v) println() println("eigenvalue") disp_eigenvalue(a) println() println("eigenvector") disp_eigenvector(v) } // ヤコビ法 def jacobi(a, v) { for (k in 1..200) { // 最大値を探す p = 0 q = 0 max_val = 0.0 for (i in 0..(N - 1)) { for (j in (i + 1)..N) { if (max_val < Math.abs(a[i][j])) { max_val = Math.abs(a[i][j]) p = i q = j } } } // θ を求める t = 0.0 if (Math.abs(a[p][p] - a[q][q]) < 0.00000000001) { // a_{pp} = a_{qq} のとき、回転角tをπ/4にする t = Math.PI / 4.0 if (a[p][p] < 0) t = -t } else { // a_{pp} ≠ a_{qq} のとき t = Math.atan(2.0 * a[p][q] / (a[p][p] - a[q][q])) / 2.0 } // θ を使って 行列 U を作成し、A = U^t × A × U c = Math.cos(t) s = Math.sin(t) // U^t × A t1 = 0.0 t2 = 0.0 for (i in 0..N) { t1 = a[p][i] * c + a[q][i] * s t2 = -a[p][i] * s + a[q][i] * c a[p][i] = t1 a[q][i] = t2 // 固有ベクトル t1 = v[p][i] * c + v[q][i] * s t2 = -v[p][i] * s + v[q][i] * c v[p][i] = t1 v[q][i] = t2 } // A × U for (i in 0..N) { t1 = a[i][p] * c + a[i][q] * s t2 = -a[i][p] * s + a[i][q] * c a[i][p] = t1 a[i][q] = t2 } // 対角要素を表示 printf("%3d\t" , k) disp_eigenvalue(a) // 収束判定 if (max_val < 0.00000000001) break } } // 対角要素を表示 def disp_eigenvalue(a) { for (i in 0..N) printf("%14.10f\t" , a[i][i]) println() } // 固有ベクトルを表示 def disp_eigenvector(matrix) { for (i in 0..N) { row = [0.0 ,0.0 ,0.0 ,0.0] for (j in 0..N) row[j] = matrix[i][j] normarize(row) disp_vector(row) } } // 1次元配列を表示 def disp_vector(row) { for (col in row) { printf("%14.10f\t" , col) } println() } // 正規化 (ベクトルの長さを1にする) def normarize(x) { s = 0.0 for (i in 0..N) s += x[i] * x[i] s = Math.sqrt(s) for (i in 0..N) x[i] /= s }
N = 3 // ハウスホルダー変換とQR分解で固有値を求める main() def main() { def a = [[5.0, 4.0, 1.0, 1.0], [4.0, 5.0, 1.0, 1.0], [1.0, 1.0, 4.0, 2.0], [1.0, 1.0, 2.0, 4.0]] as double[][] def d = [1.0 ,0.0 ,0.0 ,0.0] as double[] def e = [1.0 ,0.0 ,0.0 ,0.0] as double[] // ハウスホルダー変換 tridiagonalize(a, d, e) // QR分解 decomp(a, d, e) println() println("eigenvalue") disp_vector(d) println() println("eigenvector") disp_matrix(a) } // ハウスホルダー変換 def tridiagonalize(a, d, e) { def v = [0.0 ,0.0 ,0.0 ,0.0] as double[] for (k in 0..(N - 2)) { def w = [0.0 ,0.0 ,0.0 ,0.0] as double[] d[k] = a[k][k] t = 0.0 for (i in (k + 1)..N) { w[i] = a[i][k] t += w[i] * w[i] } s = Math.sqrt(t) if (w[k + 1] < 0) s = -s if (Math.abs(s) < 0.00000000001) e[k + 1] = 0.0 else { e[k + 1] = -s w[k + 1] += s s = 1 / Math.sqrt(w[k + 1] * s) for (i in (k + 1)..N) w[i] *= s for (i in (k + 1)..N) { s = 0.0 for (j in (k + 1)..N) { if (j <= i) s += a[i][j] * w[j] else s += a[j][i] * w[j] } v[i] = s } s = 0.0 for (i in (k + 1)..N) s += w[i] * v[i] s /= 2.0 for (i in (k + 1)..N) v[i] -= s * w[i] for (i in (k + 1)..N) for (j in (k + 1)..i) a[i][j] -= w[i] * v[j] + w[j] * v[i] // 次の行は固有ベクトルを求めないなら不要 for (i in (k + 1)..N) a[i][k] = w[i] } } d[N - 1] = a[N - 1][N - 1] d[N] = a[N][N] e[0] = 0.0 e[N] = a[N][N - 1] // 次の行は固有ベクトルを求めないなら不要 for (k in N..0) { w = [0.0 ,0.0 ,0.0 ,0.0] if (k < N - 1) { for (i in (k + 1)..N) w[i] = a[i][k] for (i in (k + 1)..N) { s = 0.0 for (j in (k + 1)..N) s += a[i][j] * w[j] v[i] = s } for (i in (k + 1)..N) for (j in (k + 1)..N) a[i][j] -= v[i] * w[j] } for (i in 0..N) a[i][k] = 0.0 a[k][k] = 1.0 } } // QR分解 def decomp(a, d, e) { e[0] = 1.0 h = N - 1 while (Math.abs(e[h]) < 0.00000000001) h-- while (h > 0) { e[0] = 0.0 l = h - 1 while (Math.abs(e[l]) >= 0.00000000001) l-- for (j in 1..100) { w = (d[h - 1] - d[h]) / 2.0 s = Math.sqrt(w * w + e[h] * e[h]) if (w < 0.0) s = -s x = d[l] - d[h] + e[h] * e[h] / (w + s) y = e[l + 1] z = 0.0 t = 0.0 u = 0.0 for (k in l..(h - 1)) { if (Math.abs(x) >= Math.abs(y)) { t = -y / x u = 1 / Math.sqrt(t * t + 1.0) s = t * u } else { t = -x / y s = 1 / Math.sqrt(t * t + 1.0) if (t < 0) s = -s u = t * s } w = d[k] - d[k + 1] t = (w * s + 2 * u * e[k + 1]) * s d[k ] = d[k ] - t d[k + 1] = d[k + 1] + t e[k ] = u * e[k] - s * z e[k + 1] = e[k + 1] * (u * u - s * s) + w * s * u // 次の行は固有ベクトルを求めないなら不要 for (i in 0..N) { x = a[k ][i] y = a[k + 1][i] a[k ][i] = u * x - s * y a[k + 1][i] = s * x + u * y } if (k < N - 1) { x = e[k + 1] y = -s * e[k + 2] z = y e[k + 2] = u * e[k + 2] } } printf("%3d\t" , j) disp_vector(d) // 収束判定 if (Math.abs(e[h]) < 0.00000000001) break } e[0] = 1.0 while (Math.abs(e[h]) < 0.00000000001) h-- } // 次の行は固有ベクトルを求めないなら不要 for (k in 0..(N - 1)) { l = k for (i in (k + 1)..N) if (d[i] > d[l]) l = i t = d[k] d[k] = d[l] d[l] = t for (i in 0..N) { t = a[k][i] a[k][i] = a[l][i] a[l][i] = t } } } // 1次元配列を表示 def disp_vector(row) { for (col in row) { printf("%14.10f\t" , col) } println() } // 2次元配列を表示 def disp_matrix(matrix) { for (row in matrix) { for (col in row) { printf("%14.10f\t" , col) } println() } }