package main import "fmt" import "math" func main() { fmt.Printf("%d\n", 3 + 5) fmt.Printf("%d\n", 3 - 5) fmt.Printf("%d\n", 3 * 5) fmt.Printf("%f\n", math.Pow(3, 5)) fmt.Printf("%d\n", 5 / 3) fmt.Printf("%f\n", 5.0 / 3) fmt.Printf("%f\n", 5 / 3.0) fmt.Printf("%d\n", 5 % 3) fmt.Printf("%3d\n", 3 * 5) fmt.Printf("%23.20f\n", 5 / 3.0) }
package main import "fmt" func main() { var i int = 3 * 5 fmt.Printf("3 * 5 = %d\n", i) }
package main import "fmt" import "math" func main() { const a = 0 const b = 1 // 台形則で積分 var n int = 2 for j := 1; j <= 10; j++ { var h float64 = (b - a) / float64(n) var s float64 = 0 var x float64 = a for i := 1; i <= n - 1; i++ { x += h s += f(x) } s = h * ((f(a) + f(b)) / 2 + s) n *= 2 // 結果を π と比較 fmt.Printf("%2d : %13.10f, %13.10f\n", j, s, s - math.Pi) } } func f(x float64) float64 { return (4 / (1 + x * x)) }
package main import "fmt" import "math" func main() { const a = 0 const b = 1 // 中点則で積分 var n int = 2 for j := 1; j <= 10; j++ { var h float64 = (b - a) / float64(n) var s float64 = 0 var x float64 = a + (h / 2) for i := 1; i <= n; i++ { s += f(x) x += h } s *= h n *= 2 // 結果を π と比較 fmt.Printf("%2d : %13.10f, %13.10f\n", j, s, s - math.Pi) } } func f(x float64) float64 { return (4 / (1 + x * x)) }
package main import "fmt" import "math" func main() { const a = 0 const b = 1 // Simpson則で積分 var n int = 2 for j := 1; j <= 5; j++ { var h float64 = (b - a) / float64(n) var s2 float64 = 0 var s4 float64 = 0 var x float64 = a + h for i := 1; i <= n / 2; i++ { s4 += f(x) x += h s2 += f(x) x += h } s2 = (s2 - f(b)) * 2 + f(a) + f(b) s4 *= 4 var s float64 = (s2 + s4) * h / 3 n *= 2 // 結果を π と比較 fmt.Printf("%2d : %13.10f, %13.10f\n", j, s, s - math.Pi) } } func f(x float64) float64 { return (4 / (1 + x * x)) }
package main import "fmt" import "math" func main() { const a = 0 const b = 1 var t[7][7] float64 // 台形則で積分 var n int = 2 for i := 1; i <= 6; i++ { var h float64 = (b - a) / float64(n) var s float64 = 0 var x float64 = a for j := 1; j <= n - 1; j++ { x += h s += f(x) } // 結果を保存 t[i][1] = h * ((f(a) + f(b)) / 2 + s) n *= 2 } // Richardsonの補外法 n = 4 for j := 2; j <= 6; j++ { for i := j; i <= 6; i++ { t[i][j] = t[i][j - 1] + (t[i][j - 1] - t[i - 1][j - 1]) / float64(n - 1) if i == j { // 結果を π と比較 fmt.Printf("%2d : %13.10f, %13.10f\n", j, t[i][j], t[i][j] - math.Pi) } } n *= 4 } } func f(x float64) float64 { return (4 / (1 + x * x)) }
package main import "fmt" import "math" // データ点の数 const N = 7 func main() { var x [N]float64 var y [N]float64 // 1.5刻みで -4.5〜4.5 まで, 7点だけ値をセット for i := 0; i < N; i++ { var d float64 = float64(i) * 1.5 - 4.5 x[i] = d y[i] = f(d) } // 0.5刻みで 与えられていない値を補間 for i := 0; i <= 18; i++ { var d float64 = float64(i) * 0.5 - 4.5 var d1 float64 = f(d) var d2 float64 = lagrange(d, x[:], y[:]) // 元の関数と比較 fmt.Printf("%5.2f\t%8.5f\t%8.5f\t%8.5f\n", d, d1, d2, d1 - d2) } } // 元の関数 func f(x float64) float64 { return x - math.Pow(x,3) / (3 * 2) + math.Pow(x,5) / (5 * 4 * 3 * 2) } // Lagrange (ラグランジュ) 補間 func lagrange(d float64, x []float64, y []float64) float64 { var sum float64 = 0 for i := 0; i < N; i++ { var prod float64 = y[i] for j := 0; j < N; j++ { if (j != i) { prod *= (d - x[j]) / (x[i] - x[j]) } } sum += prod } return sum }
package main import "fmt" import "math" // データ点の数 const N = 7 func main() { var x [N]float64 var y [N]float64 // 1.5刻みで -4.5〜4.5 まで, 7点だけ値をセット for i := 0; i < N; i++ { var d float64 = float64(i) * 1.5 - 4.5 x[i] = d y[i] = f(d) } // 0.5刻みで 与えられていない値を補間 for i := 0; i <= 18; i++ { var d float64 = float64(i) * 0.5 - 4.5 var d1 float64 = f(d) var d2 float64 = neville(d, x[:], y[:]) // 元の関数と比較 fmt.Printf("%5.2f\t%8.5f\t%8.5f\t%8.5f\n", d, d1, d2, d1 - d2) } } // 元の関数 func f(x float64) float64 { return x - math.Pow(x,3) / (3 * 2) + math.Pow(x,5) / (5 * 4 * 3 * 2) } // Neville (ネヴィル) 補間 func neville(d float64, x []float64, y []float64) float64 { var w[N][N] float64 for i := 0; i < N; i++ { w[0][i] = y[i] } for j := 1; j < N; j++ { for i := 0; i < N - j; i++ { w[j][i] = w[j-1][i+1] + (w[j-1][i+1] - w[j-1][i]) * (d - x[i+j]) / (x[i+j] - x[i]) } } return w[N-1][0] }
package main import "fmt" import "math" // データ点の数 const N = 7 func main() { var x [N]float64 var y [N]float64 // 1.5刻みで -4.5〜4.5 まで, 7点だけ値をセット for i := 0; i < N; i++ { var d float64 = float64(i) * 1.5 - 4.5 x[i] = d y[i] = f(d) } // 差分商の表を作る var d[N][N] float64 for j := 0; j < N; j++ { d[0][j] = y[j] } for i := 1; i < N; i++ { for j := 0; j < N - i; j++ { d[i][j] = (d[i-1][j+1] - d[i-1][j]) / (x[j+i] - x[j]) } } // n階差分商 var a [N]float64 for j := 0; j < N; j++ { a[j] = d[j][0] } // 0.5刻みで 与えられていない値を補間 for i := 0; i <= 18; i++ { var d float64 = float64(i) * 0.5 - 4.5 var d1 float64 = f(d) var d2 float64 = newton(d, x[:], a[:]) // 元の関数と比較 fmt.Printf("%5.2f\t%8.5f\t%8.5f\t%8.5f\n", d, d1, d2, d1 - d2) } } // 元の関数 func f(x float64) float64 { return x - math.Pow(x,3) / (3 * 2) + math.Pow(x,5) / (5 * 4 * 3 * 2) } // Newton (ニュートン) 補間 func newton(d float64, x []float64, a []float64) float64 { var sum float64 = a[0] for i := 1; i < N; i++ { var prod float64 = a[i] for j := 0; j < i; j++ { prod *= (d - x[j]) } sum += prod } return sum }
package main import "fmt" import "math" // データ点の数 const N = 7 const Nx2 = 14 func main() { var x [N]float64 var y [N]float64 var yd [N]float64 // 1.5刻みで -4.5〜4.5 まで, 7点だけ値をセット for i := 0; i < N; i++ { var d float64 = float64(i) * 1.5 - 4.5 x[i] = d y[i] = f(d) yd[i] = fd(d) } // 差分商の表を作る var z[Nx2] float64 var d[Nx2][Nx2] float64 for i := 0; i < Nx2; i++ { j := i / 2 z[i] = x[j] d[0][i] = y[j] } for i := 1; i < Nx2; i++ { for j := 0; j < Nx2 - i; j++ { if (i == 1 && j % 2 == 0) { d[i][j] = yd[j / 2] } else { d[i][j] = (d[i-1][j+1] - d[i-1][j]) / (z[j+i] - z[j]) } } } // n階差分商 var a[Nx2] float64 for j := 0; j < Nx2; j++ { a[j] = d[j][0] } // 0.5刻みで 与えられていない値を補間 for i := 0; i <= 18; i++ { var d float64 = float64(i) * 0.5 - 4.5 var d1 float64 = f(d) var d2 float64 = hermite(d, z[:], a[:]) // 元の関数と比較 fmt.Printf("%5.2f\t%8.5f\t%8.5f\t%8.5f\n", d, d1, d2, d1 - d2) } } // 元の関数 func f(x float64) float64 { return x - math.Pow(x,3) / (3 * 2) + math.Pow(x,5) / (5 * 4 * 3 * 2) } // 導関数 func fd(x float64) float64 { return 1 - math.Pow(x,2) / 2 + math.Pow(x,4) / (4 * 3 * 2) } // Hermite (エルミート) 補間 func hermite(d float64, z []float64, a []float64) float64 { var sum float64 = a[0] for i := 1; i < Nx2; i++ { var prod float64 = a[i] for j := 0; j < i; j++ { prod *= (d - z[j]) } sum += prod } return sum }
package main import "fmt" import "math" // データ点の数 const N = 7 func main() { var x [N]float64 var y [N]float64 // 1.5刻みで -4.5〜4.5 まで, 7点だけ値をセット for i := 0; i < N; i++ { var d float64 = float64(i) * 1.5 - 4.5 x[i] = d y[i] = f(d) } // 3項方程式の係数の表を作る var a [N]float64 var b [N]float64 var c [N]float64 var d [N]float64 for i := 1; i < N - 1; i++ { a[i] = x[i] - x[i-1] b[i] = 2.0 * (x[i+1] - x[i-1]) c[i] = x[i+1] - x[i] d[i] = 6.0 * ((y[i+1] - y[i]) / (x[i+1] - x[i]) - (y[i] - y[i-1]) / (x[i] - x[i-1])) } // 3項方程式を解く (ト−マス法) var g [N]float64 var s [N]float64 g[1] = b[1] s[1] = d[1] for i := 2; i < N - 1; i++ { g[i] = b[i] - a[i] * c[i-1] / g[i-1] s[i] = d[i] - a[i] * s[i-1] / g[i-1] } var z [N]float64 z[0] = 0 z[N-1] = 0 z[N-2] = s[N-2] / g[N-2] for i := N - 3; i >= 1; i-- { z[i] = (s[i] - c[i] * z[i+1]) / g[i] } // 0.5刻みで 与えられていない値を補間 for i := 0; i <= 18; i++ { var d1 float64 = float64(i) * 0.5 - 4.5 var d2 float64 = f(d1) var d3 float64 = spline(d1, x[:], y[:], z[:]) // 元の関数と比較 fmt.Printf("%5.2f\t%8.5f\t%8.5f\t%8.5f\n", d1, d2, d3, d2 - d3) } } // 元の関数 func f(x float64) float64 { return x - math.Pow(x,3) / (3 * 2) + math.Pow(x,5) / (5 * 4 * 3 * 2) } // Spline (スプライン) 補間 func spline(d float64, x []float64, y []float64, z []float64) float64 { // 補間関数値がどの区間にあるか k := -1 for i := 1; i < N; i++ { if d <= x[i] { k = i - 1 break } } if k < 0 { k = N - 1 } var d1 float64 = x[k+1] - d var d2 float64 = d - x[k] var d3 float64 = x[k+1] - x[k] return (z[k] * math.Pow(d1,3) + z[k+1] * math.Pow(d2,3)) / (6.0 * d3) + (y[k] / d3 - z[k] * d3 / 6.0) * d1 + (y[k+1] / d3 - z[k+1] * d3 / 6.0) * d2 }
package main import "fmt" import "math" // 重力加速度 const g float64 = -9.8 // 空気抵抗係数 const k float64 = -0.01 // 時間間隔(秒) const h float64 = 0.01 func main() { // 角度 var degree float64 = 45 var radian float64 = degree * math.Pi / 180.0 // 初速 250 km/h -> 秒速に変換 var v float64 = 250 * 1000 / 3600 // 水平方向の速度 var vx[2] float64 vx[0] = v * math.Cos(radian) // 鉛直方向の速度 var vy[2] float64 vy[0] = v * math.Sin(radian) // 経過秒数 var t float64 = 0.0 // 位置 var x float64 = 0.0 var y float64 = 0.0 // Euler法 for i := 1; y >= 0.0; i++ { // 経過秒数 t = float64(i) * h // 位置 x += h * vx[0] y += h * vy[0] fmt.Printf("%4.2f\t%8.5f\t%9.5f\t%9.5f\t%8.5f\n", t, vx[0], vy[0], x, y) // 速度 vx[1] = vx[0] + h * fx(vx[0], vy[0]) vy[1] = vy[0] + h * fy(vx[0], vy[0]) vx[0] = vx[1] vy[0] = vy[1] } } // 空気抵抗による水平方向の減速分 func fx(vx float64, vy float64) float64 { return k * math.Sqrt(vx * vx + vy * vy) * vx } // 重力と空気抵抗による鉛直方向の減速分 func fy(vx float64, vy float64) float64 { return g + (k * math.Sqrt(vx * vx + vy * vy) * vy) }
package main import "fmt" import "math" // 重力加速度 const g float64 = -9.8 // 空気抵抗係数 const k float64 = -0.01 // 時間間隔(秒) const h float64 = 0.01 func main() { // 角度 var degree float64 = 45 var radian float64 = degree * math.Pi / 180.0 // 初速 250 km/h -> 秒速に変換 var v float64 = 250 * 1000 / 3600 // 水平方向の速度 var vx[3] float64 vx[0] = v * math.Cos(radian) // 鉛直方向の速度 var vy[3] float64 vy[0] = v * math.Sin(radian) // 経過秒数 var t float64 = 0.0 // 位置 var x[3] float64 x[0] = 0.0 var y[3] float64 y[0] = 0.0 // Heun法 for i := 1; y[0] >= 0.0; i++ { // 経過秒数 t = float64(i) * h // 位置・速度 x[1] = x[0] + h * vx[0] y[1] = y[0] + h * vy[0] vx[1] = vx[0] + h * fx(vx[0], vy[0]) vy[1] = vy[0] + h * fy(vx[0], vy[0]) x[2] = x[0] + h * ( vx[0] + vx[1] ) / 2 y[2] = y[0] + h * ( vy[0] + vy[1] ) / 2 vx[2] = vx[0] + h * (fx(vx[0], vy[0]) + fx(vx[1], vy[1])) / 2 vy[2] = vy[0] + h * (fy(vx[0], vy[0]) + fy(vx[1], vy[1])) / 2 x[0] = x[2] y[0] = y[2] vx[0] = vx[2] vy[0] = vy[2] fmt.Printf("%4.2f\t%8.5f\t%9.5f\t%9.5f\t%9.5f\n", t, vx[0], vy[0], x[0], y[0]) } } // 空気抵抗による水平方向の減速分 func fx(vx float64, vy float64) float64 { return k * math.Sqrt(vx * vx + vy * vy) * vx } // 重力と空気抵抗による鉛直方向の減速分 func fy(vx float64, vy float64) float64 { return g + (k * math.Sqrt(vx * vx + vy * vy) * vy) }
package main import "fmt" import "math" // 重力加速度 const g float64 = -9.8 // 空気抵抗係数 const k float64 = -0.01 // 時間間隔(秒) const h float64 = 0.01 func main() { // 角度 var degree float64 = 45 var radian float64 = degree * math.Pi / 180.0 // 初速 250 km/h -> 秒速に変換 var v float64 = 250 * 1000 / 3600 // 水平方向の速度 var vx[2] float64 vx[0] = v * math.Cos(radian) // 鉛直方向の速度 var vy[2] float64 vy[0] = v * math.Sin(radian) // 経過秒数 var t float64 = 0.0 // 位置 var x[2] float64 x[0] = 0.0 var y[2] float64 y[0] = 0.0 // 中点法 for i := 1; y[0] >= 0.0; i++ { // 経過秒数 t = float64(i) * h // 位置・速度 vx[1] = h * fx(vx[0], vy[0]) vy[1] = h * fy(vx[0], vy[0]) var wx float64 = vx[0] + vx[1] / 2 var wy float64 = vy[0] + vy[1] / 2 vx[0] = vx[0] + h * fx(wx, wy) vy[0] = vy[0] + h * fy(wx, wy) x[0] = x[0] + h * wx y[0] = y[0] + h * wy fmt.Printf("%4.2f\t%8.5f\t%9.5f\t%9.5f\t%9.5f\n", t, vx[0], vy[0], x[0], y[0]) } } // 空気抵抗による水平方向の減速分 func fx(vx float64, vy float64) float64 { return k * math.Sqrt(vx * vx + vy * vy) * vx } // 重力と空気抵抗による鉛直方向の減速分 func fy(vx float64, vy float64) float64 { return g + (k * math.Sqrt(vx * vx + vy * vy) * vy) }
package main import "fmt" import "math" // 重力加速度 const g float64 = -9.8 // 空気抵抗係数 const k float64 = -0.01 // 時間間隔(秒) const h float64 = 0.01 func main() { // 角度 var degree float64 = 45 var radian float64 = degree * math.Pi / 180.0 // 初速 250 km/h -> 秒速に変換 var v float64 = 250 * 1000 / 3600 // 水平方向の速度 var vx[5] float64 vx[0] = v * math.Cos(radian) // 鉛直方向の速度 var vy[5] float64 vy[0] = v * math.Sin(radian) // 経過秒数 var t float64 = 0.0 // 位置 var x[5] float64 x[0] = 0.0 var y[5] float64 y[0] = 0.0 // Runge-Kutta法 for i := 1; y[0] >= 0.0; i++ { // 経過秒数 t = float64(i) * h // 位置・速度 x[1] = h * vx[0] y[1] = h * vy[0] vx[1] = h * fx(vx[0], vy[0]) vy[1] = h * fy(vx[0], vy[0]) var wx float64 = vx[0] + vx[1] / 2 var wy float64 = vy[0] + vy[1] / 2 x[2] = h * wx y[2] = h * wy vx[2] = h * fx(wx, wy) vy[2] = h * fy(wx, wy) wx = vx[0] + vx[2] / 2 wy = vy[0] + vy[2] / 2 x[3] = h * wx y[3] = h * wy vx[3] = h * fx(wx, wy) vy[3] = h * fy(wx, wy) wx = vx[0] + vx[3] wy = vy[0] + vy[3] x[4] = h * wx y[4] = h * wy vx[4] = h * fx(wx, wy) vy[4] = h * fy(wx, wy) x[0] += ( x[1] + x[2] * 2 + x[3] * 2 + x[4]) / 6 y[0] += ( y[1] + y[2] * 2 + y[3] * 2 + y[4]) / 6 vx[0] += (vx[1] + vx[2] * 2 + vx[3] * 2 + vx[4]) / 6 vy[0] += (vy[1] + vy[2] * 2 + vy[3] * 2 + vy[4]) / 6 fmt.Printf("%4.2f\t%8.5f\t%9.5f\t%9.5f\t%9.5f\n", t, vx[0], vy[0], x[0], y[0]) } } // 空気抵抗による水平方向の減速分 func fx(vx float64, vy float64) float64 { return k * math.Sqrt(vx * vx + vy * vy) * vx } // 重力と空気抵抗による鉛直方向の減速分 func fy(vx float64, vy float64) float64 { return g + (k * math.Sqrt(vx * vx + vy * vy) * vy) }
package main import "fmt" import "math" // 重力加速度 const g float64 = -9.8 // 空気抵抗係数 const k float64 = -0.01 // 時間間隔(秒) const h float64 = 0.01 func main() { // 角度 var degree float64 = 45 var radian float64 = degree * math.Pi / 180.0 // 初速 250 km/h -> 秒速に変換 var v float64 = 250 * 1000 / 3600 // 水平方向の速度 var vx[5] float64 vx[0] = v * math.Cos(radian) // 鉛直方向の速度 var vy[5] float64 vy[0] = v * math.Sin(radian) // 経過秒数 var t float64 = 0.0 // 位置 var x[5] float64 x[0] = 0.0 var y[5] float64 y[0] = 0.0 // Runge-Kutta-Gill法 for i := 1; y[0] >= 0.0; i++ { // 経過秒数 t = float64(i) * h // 位置・速度 x[1] = h * vx[0] y[1] = h * vy[0] vx[1] = h * fx(vx[0], vy[0]) vy[1] = h * fy(vx[0], vy[0]) var wx float64 = vx[0] + vx[1] / 2 var wy float64 = vy[0] + vy[1] / 2 x[2] = h * wx y[2] = h * wy vx[2] = h * fx(wx, wy) vy[2] = h * fy(wx, wy) wx = vx[0] + vx[1] * ((math.Sqrt(2.0) - 1) / 2) + vx[2] * (1 - 1 / math.Sqrt(2.0)) wy = vy[0] + vy[1] * ((math.Sqrt(2.0) - 1) / 2) + vy[2] * (1 - 1 / math.Sqrt(2.0)) x[3] = h * wx y[3] = h * wy vx[3] = h * fx(wx, wy) vy[3] = h * fy(wx, wy) wx = vx[0] - vx[2] / math.Sqrt(2.0) + vx[3] * (1 + 1 / math.Sqrt(2.0)) wy = vy[0] - vy[2] / math.Sqrt(2.0) + vy[3] * (1 + 1 / math.Sqrt(2.0)) x[4] = h * wx y[4] = h * wy vx[4] = h * fx(wx, wy) vy[4] = h * fy(wx, wy) x[0] += ( x[1] + x[2] * (2 - math.Sqrt(2.0)) + x[3] * (2 + math.Sqrt(2.0)) + x[4]) / 6 y[0] += ( y[1] + y[2] * (2 - math.Sqrt(2.0)) + y[3] * (2 + math.Sqrt(2.0)) + y[4]) / 6 vx[0] += (vx[1] + vx[2] * (2 - math.Sqrt(2.0)) + vx[3] * (2 + math.Sqrt(2.0)) + vx[4]) / 6 vy[0] += (vy[1] + vy[2] * (2 - math.Sqrt(2.0)) + vy[3] * (2 + math.Sqrt(2.0)) + vy[4]) / 6 fmt.Printf("%4.2f\t%8.5f\t%9.5f\t%9.5f\t%9.5f\n", t, vx[0], vy[0], x[0], y[0]) } } // 空気抵抗による水平方向の減速分 func fx(vx float64, vy float64) float64 { return k * math.Sqrt(vx * vx + vy * vy) * vx } // 重力と空気抵抗による鉛直方向の減速分 func fy(vx float64, vy float64) float64 { return g + (k * math.Sqrt(vx * vx + vy * vy) * vy) }
package main import "fmt" import "math" func main() { var a float64 = 1 var b float64 = 2 fmt.Printf("%12.10f\n", bisection(a, b)) } func bisection(a float64, b float64) float64 { var c float64 for { // 区間 (a, b) の中点 c = (a + b) / 2 c = (a + b) / 2 fmt.Printf("%12.10f\t%13.10f\n", c, c - math.Sqrt(2)) var fc float64 = f(c) if math.Fabs(fc) < 0.0000000001 { break } if fc < 0 { // f(c) < 0 であれば, 解は区間 (c, b) の中に存在 a = c } else { // f(c) > 0 であれば, 解は区間 (a, c) の中に存在 b = c } } return c } func f(x float64) float64 { return x * x - 2 }
package main import "fmt" import "math" func main() { var a float64 = 1 var b float64 = 2 fmt.Printf("%12.10f\n", falseposition(a, b)) } func falseposition(a float64, b float64) float64 { var c float64 for { // 点 (a,f(a)) と 点 (b,f(b)) を結ぶ直線と x軸の交点 c = (a * f(b) - b * f(a)) / (f(b) - f(a)) fmt.Printf("%12.10f\t%13.10f\n", c, c - math.Sqrt(2)) var fc float64 = f(c) if math.Fabs(fc) < 0.0000000001 { break } if fc < 0 { // f(c) < 0 であれば, 解は区間 (c, b) の中に存在 a = c } else { // f(c) > 0 であれば, 解は区間 (a, c) の中に存在 b = c } } return c } func f(x float64) float64 { return x * x - 2 }
package main import "fmt" import "math" func main() { var x float64 = 1 fmt.Printf("%12.10f\n", iterative(x)) } func iterative(x0 float64) float64 { var x1 float64 for { x1 = g(x0) fmt.Printf("%12.10f\t%13.10f\n", x1, x1 - math.Sqrt(2)) if math.Fabs(x1 - x0) < 0.0000000001 { break } x0 = x1 } return x1 } func g(x float64) float64 { return (x / 2) + (1 / x) }
package main import "fmt" import "math" func main() { var x float64 = 2 fmt.Printf("%12.10f\n", newton(x)) } func newton(x0 float64) float64 { var x1 float64 for { x1 = x0 - (f0(x0) / f1(x0)) fmt.Printf("%12.10f\t%13.10f\n", x1, x1 - math.Sqrt(2)) if math.Fabs(x1 - x0) < 0.0000000001 { break } x0 = x1 } return x1 } func f0(x float64) float64 { return x * x - 2 } func f1(x float64) float64 { return 2 * x }
package main import "fmt" import "math" func main() { var x float64 = 2 fmt.Printf("%12.10f\n", bailey(x)) } func bailey(x0 float64) float64 { var x1 float64 for { x1 = x0 - (f0(x0) / (f1(x0) - (f0(x0) * f2(x0) / (2 * f1(x0))))) fmt.Printf("%12.10f\t%13.10f\n", x1, x1 - math.Sqrt(2)) if math.Fabs(x1 - x0) < 0.0000000001 { break } x0 = x1 } return x1 } func f0(x float64) float64 { return x * x - 2 } func f1(x float64) float64 { return 2 * x } func f2(x float64) float64 { return 2 }
package main import "fmt" import "math" func main() { var x0 float64 = 1 var x1 float64 = 2 fmt.Printf("%12.10f\n", secant(x0, x1)) } func secant(x0 float64, x1 float64) float64 { var x2 float64 for { x2 = x1 - f(x1) * (x1 - x0) / (f(x1) - f(x0)) fmt.Printf("%12.10f\t%13.10f\n", x2, x2 - math.Sqrt(2)) if math.Fabs(x2 - x1) < 0.0000000001 { break } x0 = x1 x1 = x2 } return x2 } func f(x float64) float64 { return x * x - 2 }
package main import "fmt" import "math" const N = 4 func main() { var a [N][N]float64 = [N][N]float64{{9,2,1,1},{2,8,-2,1},{-1,-2,7,-2},{1,-1,-2,6}} var b = []float64{20,16,8,17} var c = []float64{0,0,0,0} // ヤコビの反復法 jacobi(a,b,c) fmt.Println("X") disp_vector(c) } // ヤコビの反復法 func jacobi(a[N][N]float64, b[]float64, x0[]float64) { for { var x1 = []float64{0,0,0,0} var finish = true for i := 0; i < N; i++ { x1[i] = b[i] for j := 0; j < N; j++ { if (j != i) { x1[i] -= a[i][j] * x0[j] } } x1[i] /= a[i][i] if (math.Fabs(x1[i] - x0[i]) > 0.0000000001) { finish = false } } for i := 0; i < N; i++ { x0[i] = x1[i] } if (finish) { return } disp_vector(x0) } } // 1次元配列を表示 func disp_vector(row[]float64) { for _, col := range row { fmt.Printf("%14.10f\t", col) } fmt.Println("") }
package main import "fmt" import "math" const N = 4 func main() { var a [N][N]float64 = [N][N]float64{{9,2,1,1},{2,8,-2,1},{-1,-2,7,-2},{1,-1,-2,6}} var b = []float64{20,16,8,17} var c = []float64{0,0,0,0} // ガウス・ザイデル法 gauss(a,b,c) fmt.Println("X") disp_vector(c) } // ガウス・ザイデル法 func gauss(a[N][N]float64, b[]float64, x0[]float64) { for { var finish = true for i := 0; i < N; i++ { var x1 float64 = b[i] for j := 0; j < N; j++ { if (j != i) { x1 -= a[i][j] * x0[j] } } x1 /= a[i][i] if (math.Fabs(x1 - x0[i]) > 0.0000000001) { finish = false } x0[i] = x1 } if (finish) { return } disp_vector(x0) } } // 1次元配列を表示 func disp_vector(row[]float64) { for _, col := range row { fmt.Printf("%14.10f\t", col) } fmt.Println("") }
package main import "fmt" import "math" const N = 4 // ガウスの消去法 func main() { var a [N][N]float64 = [N][N]float64{{-1,-2,7,-2},{1,-1,-2,6},{9,2,1,1},{2,8,-2,1}} var b []float64 = []float64{8,17,20,16} // ピボット選択 pivoting(&a,b) fmt.Println("pivoting") fmt.Println("A") disp_matrix(a); fmt.Println("B") disp_vector(b); fmt.Println() // 前進消去 forward_elimination(&a,b) fmt.Println("forward elimination") fmt.Println("A") disp_matrix(a); fmt.Println("B") disp_vector(b); fmt.Println() // 後退代入 backward_substitution(&a,b) fmt.Println("X") disp_vector(b); } // 前進消去 func forward_elimination(a *[N][N]float64, b []float64) { for pivot := 0; pivot < N - 1; pivot++ { for row := pivot + 1; row < N; row++ { var s = a[row][pivot] / a[pivot][pivot] for col := pivot; col < N; col++ { a[row][col] -= a[pivot][col] * s } b[row] -= b[pivot] * s } } } // 後退代入 func backward_substitution(a *[N][N]float64, b []float64) { for row := N - 1; row >= 0; row-- { for col := N - 1; col > row; col-- { b[row] -= a[row][col] * b[col] } b[row] /= a[row][row] } } // ピボット選択 func pivoting(a *[N][N]float64, b []float64) { for pivot := 0; pivot < N; pivot++ { // 各列で 一番値が大きい行を 探す var max_row = pivot var max_val float64 = 0.0 for row := pivot; row < N; row++ { if math.Fabs(a[row][pivot]) > max_val { // 一番値が大きい行 max_val = math.Fabs(a[row][pivot]) max_row = row } } // 一番値が大きい行と入れ替え if (max_row != pivot) { var tmp float64 = 0.0 for col := 0; col < N; col++ { 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次元配列を表示 func disp_vector(row[]float64) { for _, col := range row { fmt.Printf("%14.10f\t", col) } fmt.Println("") } // 2次元配列を表示 func disp_matrix(matrix[N][N]float64) { for _, row := range matrix { for _, col := range row { fmt.Printf("%14.10f\t", col) } fmt.Println("") } }
package main import "fmt" import "math" const N = 4 // ガウス・ジョルダン法 func main() { var a [N][N]float64 = [N][N]float64{{-1,-2,7,-2},{1,-1,-2,6},{9,2,1,1},{2,8,-2,1}} var b []float64 = []float64{8,17,20,16} // ピボット選択 pivoting(&a,b) fmt.Println("pivoting") fmt.Println("A") disp_matrix(a); fmt.Println("B") disp_vector(b); fmt.Println() // 前進消去 forward_elimination(&a,b) fmt.Println("forward elimination") fmt.Println("A") disp_matrix(a); fmt.Println("B") disp_vector(b); fmt.Println() // 後退代入 backward_substitution(&a,b) fmt.Println("X") disp_vector(b); } // 前進消去 func forward_elimination(a *[N][N]float64, b []float64) { for pivot := 0; pivot < N; pivot++ { for row := 0; row < N; row++ { if (row == pivot) { continue } var s = a[row][pivot] / a[pivot][pivot] for col := pivot; col < N; col++ { a[row][col] -= a[pivot][col] * s } b[row] -= b[pivot] * s } } } // 後退代入 func backward_substitution(a *[N][N]float64, b []float64) { for pivot := 0; pivot < N; pivot++ { b[pivot] /= a[pivot][pivot] } } // ピボット選択 func pivoting(a *[N][N]float64, b []float64) { for pivot := 0; pivot < N; pivot++ { // 各列で 一番値が大きい行を 探す var max_row = pivot var max_val float64 = 0.0 for row := pivot; row < N; row++ { if math.Fabs(a[row][pivot]) > max_val { // 一番値が大きい行 max_val = math.Fabs(a[row][pivot]) max_row = row } } // 一番値が大きい行と入れ替え if (max_row != pivot) { var tmp float64 = 0.0 for col := 0; col < N; col++ { 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次元配列を表示 func disp_vector(row[]float64) { for _, col := range row { fmt.Printf("%14.10f\t", col) } fmt.Println("") } // 2次元配列を表示 func disp_matrix(matrix[N][N]float64) { for _, row := range matrix { for _, col := range row { fmt.Printf("%14.10f\t", col) } fmt.Println("") } }
package main import "fmt" import "math" const N = 4 // LU分解 func main() { var a [N][N]float64 = [N][N]float64{{-1,-2,7,-2},{1,-1,-2,6},{9,2,1,1},{2,8,-2,1}} var b []float64 = []float64{8,17,20,16} // ピボット選択 pivoting(&a,b) fmt.Println("pivoting") fmt.Println("A") disp_matrix(a) fmt.Println("B") disp_vector(b) fmt.Println() // 前進消去 forward_elimination(&a,b) fmt.Println("LU") disp_matrix(a) // Ly=b から y を求める (前進代入) var y []float64 = []float64{0,0,0,0} forward_substitution(a,b,y) fmt.Println("Y") disp_vector(y) // Ux=y から x を求める (後退代入) var x []float64 = []float64{0,0,0,0} backward_substitution(a,y,x) fmt.Println("X") disp_vector(x) } // 前進消去 func forward_elimination(a *[N][N]float64, b []float64) { for pivot := 0; pivot < N - 1; pivot++ { for row := pivot + 1; row < N; row++ { var s = a[row][pivot] / a[pivot][pivot] for col := pivot; col < N; col++ { a[row][col] -= a[pivot][col] * s // これが 上三角行列 } a[row][pivot] = s // これが 下三角行列 // b[row] -= b[pivot] * s // この値は変更しない } } } // 前進代入 func forward_substitution(a [N][N]float64, b []float64, y []float64) { for row := 0; row < N; row++ { for col := 0; col < row; col++ { b[row] -= a[row][col] * y[col] } y[row] = b[row] } } // 後退代入 func backward_substitution(a [N][N]float64, y []float64, x []float64) { for row := N - 1; row >= 0; row-- { for col := N - 1; col > row; col-- { y[row] -= a[row][col] * x[col] } x[row] = y[row] / a[row][row] } } // ピボット選択 func pivoting(a *[N][N]float64, b []float64) { for pivot := 0; pivot < N; pivot++ { // 各列で 一番値が大きい行を 探す var max_row = pivot var max_val float64 = 0.0 for row := pivot; row < N; row++ { if math.Fabs(a[row][pivot]) > max_val { // 一番値が大きい行 max_val = math.Fabs(a[row][pivot]) max_row = row } } // 一番値が大きい行と入れ替え if (max_row != pivot) { var tmp float64 = 0.0 for col := 0; col < N; col++ { 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次元配列を表示 func disp_vector(row[]float64) { for _, col := range row { fmt.Printf("%14.10f\t", col) } fmt.Println("") } // 2次元配列を表示 func disp_matrix(matrix[N][N]float64) { for _, row := range matrix { for _, col := range row { fmt.Printf("%14.10f\t", col) } fmt.Println("") } }
package main import "fmt" import "math" const N = 4 // コレスキー法 func main() { var a [N][N]float64 = [N][N]float64{{5,2,3,4},{2,10,6,7},{3,6,15,9},{4,7,9,20}} var b []float64 = []float64{34,68,96,125} fmt.Println("A") disp_matrix(a) fmt.Println("B") disp_vector(b) fmt.Println() // 前進消去 forward_elimination(&a,b) fmt.Println("LL^T") disp_matrix(a) // Ly=b から y を求める (前進代入) var y []float64 = []float64{0,0,0,0} forward_substitution(a,b,y) fmt.Println("Y") disp_vector(y) // L^Tx=y から x を求める (後退代入) var x []float64 = []float64{0,0,0,0} backward_substitution(a,y,x) fmt.Println("X") disp_vector(x) } // 前進消去 func forward_elimination(a *[N][N]float64, b []float64) { for pivot := 0; pivot < N; pivot++ { var s = 0.0 for col := 0; col < pivot; col++ { s += a[pivot][col] * a[pivot][col] } // ここで根号の中が負の値になると計算できない! a[pivot][pivot] = math.Sqrt(a[pivot][pivot] - s) for row := pivot + 1; row < N; row++ { s = a[row][pivot] / a[pivot][pivot] s = 0 for col := 0; col < pivot; col++ { s += a[row][col] * a[pivot][col] } a[row][pivot] = (a[row][pivot] - s) / a[pivot][pivot] a[pivot][row] = a[row][pivot] } } } // 前進代入 func forward_substitution(a [N][N]float64, b []float64, y []float64) { for row := 0; row < N; row++ { for col := 0; col < row; col++ { b[row] -= a[row][col] * y[col] } y[row] = b[row] / a[row][row] } } // 後退代入 func backward_substitution(a [N][N]float64, y []float64, x []float64) { for row := N - 1; row >= 0; row-- { for col := N - 1; col > row; col-- { y[row] -= a[col][row] * x[col] } x[row] = y[row] / a[row][row] } } // 1次元配列を表示 func disp_vector(row[]float64) { for _, col := range row { fmt.Printf("%14.10f\t", col) } fmt.Println("") } // 2次元配列を表示 func disp_matrix(matrix[N][N]float64) { for _, row := range matrix { for _, col := range row { fmt.Printf("%14.10f\t", col) } fmt.Println("") } }
package main import "fmt" const N = 4 // 修正コレスキー法 func main() { var a [N][N]float64 = [N][N]float64{{5,2,3,4},{2,10,6,7},{3,6,15,9},{4,7,9,20}} var b []float64 = []float64{34,68,96,125} fmt.Println("A") disp_matrix(a) fmt.Println("B") disp_vector(b) fmt.Println() // 前進消去 forward_elimination(&a,b) fmt.Println("LDL^T") disp_matrix(a) // Ly=b から y を求める (前進代入) var y []float64 = []float64{0,0,0,0} forward_substitution(a,b,y) fmt.Println("Y") disp_vector(y) // DL^Tx=y から x を求める (後退代入) var x []float64 = []float64{0,0,0,0} backward_substitution(a,y,x) fmt.Println("X") disp_vector(x) } // 前進消去 func forward_elimination(a *[N][N]float64, b []float64) { for pivot := 0; pivot < N; pivot++ { var s = 0.0 // pivot < k の場合 for col := 0; col < pivot; col++ { s = a[pivot][col] for k := 0; k < col; k++ { 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] for k := 0; k < pivot; k++ { s -= a[pivot][k] * a[pivot][k] * a[k][k] } a[pivot][pivot] = s } } // 前進代入 func forward_substitution(a [N][N]float64, b []float64, y []float64) { for row := 0; row < N; row++ { for col := 0; col < row; col++ { b[row] -= a[row][col] * y[col] } y[row] = b[row] } } // 後退代入 func backward_substitution(a [N][N]float64, y []float64, x []float64) { for row := N - 1; row >= 0; row-- { for col := N - 1; col > row; col-- { y[row] -= a[col][row] * a[row][row] * x[col] } x[row] = y[row] / a[row][row] } } // 1次元配列を表示 func disp_vector(row[]float64) { for _, col := range row { fmt.Printf("%14.10f\t", col) } fmt.Println("") } // 2次元配列を表示 func disp_matrix(matrix[N][N]float64) { for _, row := range matrix { for _, col := range row { fmt.Printf("%14.10f\t", col) } fmt.Println("") } }
package main import "fmt" import "math" const N = 4 // ベキ乗法で最大固有値を求める func main() { var a [N][N]float64 = [N][N]float64{{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}} var x = []float64 {1.0 ,0.0 ,0.0 ,0.0} // ベキ乗法 var lambda float64 = power(a, x) fmt.Println("\neigenvalue") fmt.Printf("%14.10f\n", lambda) fmt.Println("eigenvector") disp_vector(x) } // ベキ乗法 func power(a[N][N]float64, x0[]float64) float64 { var lambda float64 = 0.0 // 正規化 (ベクトル x0 の長さを1にする) normarize(x0) var e0 float64 = 0.0 for i := 0; i < N; i++ { e0 += x0[i] } for k := 1; k < 200; k++ { // 1次元配列を表示 fmt.Printf("%3d\t", k) disp_vector(x0) // 行列の積 x1 = A × x0 var x1 = []float64 {0.0 ,0.0 ,0.0 ,0.0} for i := 0; i < N; i++ { for j := 0; j < N; j++ { x1[i] += a[i][j] * x0[j] } } // 内積 var p0 float64 = 0.0 var p1 float64 = 0.0 for i := 0; i < N; i++ { p0 += x1[i] * x1[i] p1 += x1[i] * x0[i] } // 固有値 lambda = p0 / p1 // 正規化 (ベクトル x1 の長さを1にする) normarize(x1) var e1 float64 = 0.0 // 収束判定 for i := 0; i < N; i++ { e1 += x1[i] } if (math.Abs(e0 - e1) < 0.00000000001) { break } for i := 0; i < N; i++ { x0[i] = x1[i] } e0 = e1 } return lambda } // 1次元配列を表示 func disp_vector(row[]float64) { for _, col := range row { fmt.Printf("%14.10f\t", col) } fmt.Println("") } // 正規化 (ベクトルの長さを1にする) func normarize(x[]float64) { var s float64 = 0.0 for i := 0; i < N; i++ { s += x[i] * x[i] } s = math.Sqrt(s) for i := 0; i < N; i++ { x[i] /= s } }
package main import "fmt" import "math" const N = 4 func main() { var a [N][N]float64 = [N][N]float64{{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}} var x = []float64 {1.0 ,0.0 ,0.0 ,0.0} // LU分解 forward_elimination(&a) // 逆ベキ乗法 var lambda float64 = inverse(a, x) fmt.Println("\neigenvalue") fmt.Printf("%14.10f\n", lambda) fmt.Println("eigenvector") disp_vector(x) } // 逆ベキ乗法 func inverse(a[N][N]float64, x0[]float64) float64 { var lambda float64 = 0.0 // 正規化 (ベクトル x0 の長さを1にする) normarize(x0) var e0 float64 = 0.0 for i := 0; i < N; i++ { e0 += x0[i] } for k := 1; k < 200; k++ { // 1次元配列を表示 fmt.Printf("%3d\t", k) disp_vector(x0) // Ly = b から y を求める (前進代入) var b = []float64 {0.0 ,0.0 ,0.0 ,0.0} var y = []float64 {0.0 ,0.0 ,0.0 ,0.0} for i := 0; i < N; i++ { b[i] = x0[i] } forward_substitution(a,y,b) // Ux = y から x を求める (後退代入) var x1 = []float64 {0.0 ,0.0 ,0.0 ,0.0} backward_substitution(a,x1,y) // 内積 var p0 float64 = 0.0 var p1 float64 = 0.0 for i := 0; i < N; i++ { p0 += x1[i] * x1[i] p1 += x1[i] * x0[i] } // 固有値 lambda = p1 / p0 // 正規化 (ベクトル x1 の長さを1にする) normarize(x1) // 収束判定 var e1 float64 = 0.0 for i := 0; i < N; i++ { e1 += x1[i] } if (math.Abs(e0 - e1) < 0.00000000001) { break } for i := 0; i < N; i++ { x0[i] = x1[i] } e0 = e1 } return lambda } // LU分解 func forward_elimination(a *[N][N]float64) { for pivot := 0; pivot < N - 1; pivot++ { for row := pivot + 1; row < N; row++ { var s = a[row][pivot] / a[pivot][pivot] for col := pivot; col < N; col++ { a[row][col] -= a[pivot][col] * s // これが 上三角行列 } a[row][pivot] = s // これが 下三角行列 } } } // 前進代入 func forward_substitution(a [N][N]float64, y []float64, b []float64) { for row := 0; row < N; row++ { for col := 0; col < row; col++ { b[row] -= a[row][col] * y[col] } y[row] = b[row] } } // 後退代入 func backward_substitution(a [N][N]float64, x []float64, y []float64) { for row := N - 1; row >= 0; row-- { for col := N - 1; col > row; col-- { y[row] -= a[row][col] * x[col] } x[row] = y[row] / a[row][row] } } // 1次元配列を表示 func disp_vector(row[]float64) { for _, col := range row { fmt.Printf("%14.10f\t", col) } fmt.Println("") } // 正規化 (ベクトルの長さを1にする) func normarize(x[]float64) { var s float64 = 0.0 for i := 0; i < N; i++ { s += x[i] * x[i] } s = math.Sqrt(s) for i := 0; i < N; i++ { x[i] /= s } }
package main import "fmt" import "math" const N = 4 // LR分解で固有値を求める func main() { var a [N][N]float64 = [N][N]float64{{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}} var l [N][N]float64 = [N][N]float64{{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}} var u [N][N]float64 = [N][N]float64{{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}} for k := 1; k < 200; k++ { // LU分解 decomp(a, &l, &u) // 行列の積 multiply(u, l, &a) // 対角要素を表示 fmt.Printf("%3d\t", k) disp_eigenvalue(a) // 収束判定 var e float64 = 0.0 for i := 1; i < N; i++ { for j := 0; j < i; j++ { e += math.Abs(a[i][j]) } } if (e < 0.00000000001) { break } } fmt.Println("\neigenvector") disp_eigenvalue(a) } // LU分解 func decomp(a[N][N]float64, l *[N][N]float64, u *[N][N]float64) { for i := 0; i < N; i++ { for j := 0; j < N; j++ { l[i][j] = 0.0 u[i][j] = 0.0 } } l[0][0] = 1.0 for j := 0; j < N; j++ { u[0][j] = a[0][j] } for i := 1; i < N; i++ { u[i][0] = 0.0 l[0][i] = 0.0 l[i][0] = a[i][0] / u[0][0] } for i := 1; i < N; i++ { l[i][i] = 1.0 var t float64 = a[i][i] for k := 0; k <= i; k++ { t -= l[i][k] * u[k][i] } u[i][i] = t for j := i + 1; j < N; j++ { u[j][i] = 0.0 l[i][j] = 0.0 t = a[j][i] for k := 0; k <= i; k++ { t -= l[j][k] * u[k][i] } l[j][i] = t / u[i][i] t = a[i][j] for k := 0; k <= i; k++ { t -= l[i][k] * u[k][j] } u[i][j] = t } } } // 行列の積 func multiply(a[N][N]float64, b[N][N]float64, c *[N][N]float64) { for i := 0; i < N; i++ { for j := 0; j < N; j++ { var s float64 = 0.0 for k := 0; k < N; k++ { s += a[i][k] * b[k][j] } c[i][j] = s } } } // 対角要素を表示 func disp_eigenvalue(a[N][N]float64) { for i := 0; i < N; i++ { fmt.Printf("%14.10f\t", a[i][i]) } fmt.Println("") }
package main import "fmt" import "math" const N = 4 // QR分解で固有値を求める func main() { var a [N][N]float64 = [N][N]float64{{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}} var q [N][N]float64 = [N][N]float64{{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}} var r [N][N]float64 = [N][N]float64{{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}} for k := 1; k < 200; k++ { // QR分解 decomp(a, &q, &r) // 行列の積 multiply(r, q, &a) // 対角要素を表示 fmt.Printf("%3d\t", k) disp_eigenvalue(a) // 収束判定 var e float64 = 0.0 for i := 1; i < N; i++ { for j := 0; j < i; j++ { e += math.Abs(a[i][j]) } } if (e < 0.00000000001) { break } } fmt.Println("\neigenvalue") disp_eigenvalue(a) } // QR分解 func decomp(a[N][N]float64, q *[N][N]float64, r *[N][N]float64) { var x = []float64 {0.0 ,0.0 ,0.0 ,0.0} for k := 0; k < N; k++ { for i := 0; i < N; i++ { x[i] = a[i][k] } for j := 0; j < k; j++ { var t float64 = 0.0 for i := 0; i < N; i++ { t += a[i][k] * q[i][j] } r[j][k] = t r[k][j] = 0.0 for i := 0; i < N; i++ { x[i] -= t * q[i][j] } } var s float64 = 0.0 for i := 0; i < N; i++ { s += x[i] * x[i] } r[k][k] = math.Sqrt(s) for i := 0; i < N; i++ { q[i][k] = x[i] / r[k][k] } } } // 行列の積 func multiply(a[N][N]float64, b[N][N]float64, c *[N][N]float64) { for i := 0; i < N; i++ { for j := 0; j < N; j++ { var s float64 = 0.0 for k := 0; k < N; k++ { s += a[i][k] * b[k][j] } c[i][j] = s } } } // 対角要素を表示 func disp_eigenvalue(a[N][N]float64) { for i := 0; i < N; i++ { fmt.Printf("%14.10f\t", a[i][i]) } fmt.Println("") }
package main import "fmt" import "math" const N = 4 // ヤコビ法で固有値を求める func main() { var a [N][N]float64 = [N][N]float64{{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}} var v [N][N]float64 = [N][N]float64{{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}} // ヤコビ法 jacobi(&a, &v) fmt.Println("\neigenvalue") disp_eigenvalue(a) fmt.Println("\neigenvector") disp_eigenvector(v) } // ヤコビ法 func jacobi(a *[N][N]float64, v *[N][N]float64) { for k := 1; k < 100; k++ { // 最大値を探す var p int = 0 var q int = 0 var max_val float64 = 0.0 for i := 0; i < N; i++ { for j := i + 1; j < N; j++ { if (max_val < math.Abs(a[i][j])) { max_val = math.Abs(a[i][j]) p = i q = j } } } // θ を求める var t float64 = 0.0 if (math.Abs(a[p][p] - a[q][q]) < 0.00000000001) { // a_{pp} = a_{qq} のとき、回転角θをπ/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 } // θ を使って 行列 P を作成し、A = P^t × A × P var c float64 = math.Cos(t) var s float64 = math.Sin(t) // P^t × A var t1 float64 = 0.0 var t2 float64 = 0.0 for i := 0; i < N; i++ { 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 × P for i := 0; i < N; i++ { 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 } // 対角要素を表示 fmt.Printf("%3d\t", k) disp_eigenvalue(*a) // 収束判定 if (max_val < 0.00000000001) { break } } } // 対角要素を表示 func disp_eigenvalue(a[N][N]float64) { for i := 0; i < N; i++ { fmt.Printf("%14.10f\t", a[i][i]) } fmt.Println("") } // 固有ベクトルを表示 func disp_eigenvector(matrix[N][N]float64) { for i := 0; i < N; i++ { var row = []float64 {0.0 ,0.0 ,0.0 ,0.0} for j := 0; j < N; j++ { row[j] = matrix[i][j] } normarize(row) disp_vector(row) } } // 1次元配列を表示 func disp_vector(row[]float64) { for _, col := range row { fmt.Printf("%14.10f\t", col) } fmt.Println("") } // 正規化 (ベクトルの長さを1にする) func normarize(x[]float64) { var s float64 = 0.0 for i := 0; i < N; i++ { s += x[i] * x[i] } s = math.Sqrt(s) for i := 0; i < N; i++ { x[i] /= s } }
package main import "fmt" import "math" const N = 4 // ハウスホルダー変換とQR分解で固有値を求める func main() { var a [N][N]float64 = [N][N]float64{{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}} var d = []float64 {0.0 ,0.0 ,0.0 ,0.0} var e = []float64 {0.0 ,0.0 ,0.0 ,0.0} // ハウスホルダー変換 tridiagonalize(&a, d, e) // QR分解 decomp(&a, d, e) fmt.Println("\neigenvalue") disp_vector(d) fmt.Println("\neigenvector") disp_matrix(a) } // ハウスホルダー変換 func tridiagonalize(a *[N][N]float64, d []float64, e []float64) { var v = []float64 {0.0 ,0.0 ,0.0 ,0.0} for k := 0; k < N - 2; k++ { var w = []float64 {0.0 ,0.0 ,0.0 ,0.0} d[k] = a[k][k] var t float64 = 0.0 for i := k + 1; i < N; i++ { w[i] = a[i][k] t += w[i] * w[i] } var s float64 = 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 := k + 1; i < N; i++ { w[i] *= s } for i := k + 1; i < N; i++ { s = 0.0 for j := k + 1; j < N; j++ { if (j <= i) { s += a[i][j] * w[j] } else { s += a[j][i] * w[j] } } v[i] = s } s = 0.0 for i := k + 1; i < N; i++ { s += w[i] * v[i] } s /= 2.0 for i := k + 1; i < N; i++ { v[i] -= s * w[i] } for i := k + 1; i < N; i++ { for j := k + 1; j < i + 1; j++ { a[i][j] -= w[i] * v[j] + w[j] * v[i] } } // 次の行は固有ベクトルを求めないなら不要 for i := k + 1; i < N; i++ { a[i][k] = w[i] } } } d[N - 2] = a[N - 2][N - 2] d[N - 1] = a[N - 1][N - 1] e[0] = 0.0 e[N - 1] = a[N - 1][N - 2] // 次の行は固有ベクトルを求めないなら不要 for k := N - 1; k >= 0; k-- { var w = []float64 {0.0 ,0.0 ,0.0 ,0.0} if (k < N - 2) { for i := k + 1; i < N; i++ { w[i] = a[i][k] } for i := k + 1; i < N; i++ { var s float64 = 0.0 for j := k + 1; j < N; j++ { s += a[i][j] * w[j] } v[i] = s } for i := k + 1; i < N; i++ { for j := k + 1; j < N; j++ { a[i][j] -= v[i] * w[j] } } } for i := 0; i < N; i++ { a[i][k] = 0.0 } a[k][k] = 1.0 } } // QR分解 func decomp(a *[N][N]float64, d []float64, e []float64) { e[0] = 1.0 var h int = N - 1 for (math.Abs(e[h]) < 0.00000000001) { h-- } for (h > 0) { e[0] = 0.0 var l int = h - 1 for (math.Abs(e[l]) >= 0.00000000001) { l-- } for j := 1; j < 100; j++ { var w float64 = (d[h - 1] - d[h]) / 2.0 var s float64 = math.Sqrt(w * w + e[h] * e[h]) if (w < 0.0) { s = -s } var x float64 = d[l] - d[h] + e[h] * e[h] / (w + s) var y float64 = e[l + 1] var z float64 = 0.0 var t float64 = 0.0 var u float64 = 0.0 for k := l; k < h; k++ { 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 := 0; i < N; i++ { 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 - 2) { x = e[k + 1] y = -s * e[k + 2] z = y e[k + 2] = u * e[k + 2] } } fmt.Printf("%3d\t", j) disp_vector(d) // 収束判定 if (math.Abs(e[h]) < 0.00000000001) { break } } e[0] = 1.0 for (math.Abs(e[h]) < 0.00000000001) { h-- } } // 次の行は固有ベクトルを求めないなら不要 for k := 0; k < N - 1; k++ { var l int = k for i := k + 1; i < N; i++ { if (d[i] > d[l]) { l = i } } var t float64 = d[k] d[k] = d[l] d[l] = t for i := 0; i < N; i++ { t = a[k][i] a[k][i] = a[l][i] a[l][i] = t } } } // 1次元配列を表示 func disp_vector(row[]float64) { for _, col := range row { fmt.Printf("%14.10f\t", col) } fmt.Println("") } // 2次元配列を表示 func disp_matrix(matrix[N][N]float64) { for row := 0; row < N; row++ { for col := 0; col < N; col++ { fmt.Printf("%14.10f\t", matrix[row][col]) } fmt.Println() } }