puts 3 + 5 puts 3 - 5 puts 3 * 5 puts 3 ** 5 puts 5 / 3.0 puts 5.0 / 3 puts 5 / 3 puts 5 % 3 print 3 * 5, "\n" printf "%3d\n", 3 * 5 printf "%23.20f\n", 5 / 3.0
i = 3 * 5 print "3 * 5 = ", i, "\n" print "3 * 5 = #{i}\n" printf "3 * 5 = %d\n", i puts "3 * 5 = #{i}"
(1..9).each do |i| print "#{i}, " end print "\n"
(1..9).each do |i| if i % 3 == 0 print "#{i}, " end end print "\n" (1..9).each do |i| print "#{i}, " if i % 3 == 0 end print "\n" (1..9).each do |i| print "#{i}, " unless i % 3 != 0 end print "\n"
sum = 0 (1..99).each do |i| if i % 3 == 0 sum += i end end puts sum
1..9 (1..9).each { |i| print i, ", " }
(1..9).select { |n| n % 3 == 0}
(1..99).select { |n| n % 3 == 0}.inject(:+) (1..99).select { |n| n % 3 == 0}.inject("+")
# 初項:a, 公差:a で 上限:lim の数列の総和を返す関数 def sn(a, lim) n = lim / a # 項数:n = 上限:lim / 公差:a l = n * a # 末項:l = 項数:n * 公差:a (a + l) * n / 2 # 総和:sn = (初項:a + 末項:l) * 項数:n / 2 end # 3 の倍数の合計 sum = sn(3, 999) puts sum
# 10000 までの 自然数の和 # 項目数 n = 10000 n = 10000 puts (n * (n + 1) / 2)
# 10000 までの 偶数の和 # 項目数 n = 5000 n = 10000 / 2 puts (n * (n + 1))
# 10000 までの 奇数の和 # 項目数 n = 5000 n = 10000 / 2 puts (n ** 2)
# 1000 までの 自然数の2乗の和 n = 1000 puts (n * (n + 1) * (2 * n + 1) / 6)
# 100 までの 自然数の3乗の和 n = 100 puts ((n ** 2) * ((n + 1) ** 2) / 4)
# 初項 2, 公比 3, 項数 10 の等比数列の和 n = 10 a = 2 r = 3 puts ((a * ((r ** n) - 1)) / (r - 1))
a = 5 # 初項 5 d = 3 # 公差 3 n = 10 # 項数 10 p = 1 # 積 (1..n).each do |i| m = a + (d * (i - 1)) p *= m end puts p
def product(m, d, n) if (n == 0) 1 else m * product(m + d, d, n - 1) end end # 初項 5, 公差 3, 項数 10 の数列の積を表示する puts product(5, 3, 10)
# 階乗を求める関数 def Fact(n) if (n <= 1) 1 else n * Fact(n - 1) end end # 10の階乗 puts Fact(10) puts 10 * 9 * 8 * 7 * 6 * 5 * 4 * 3 * 2 * 1
# 下降階乗冪 def FallingFact(x, n) if (n <= 1) x else x * FallingFact(x - 1, n - 1) end end # 10 から 6 までの 総乗 puts FallingFact(10, 5) puts 10 * 9 * 8 * 7 * 6
# 上昇階乗冪 def RisingFact(x, n) if (n <= 1) x else x * RisingFact(x + 1, n - 1) end end # 10 から 14 までの 総乗 puts RisingFact(10, 5) puts 10 * 11 * 12 * 13 * 14
# 階乗 def Fact(n) if (n <= 1) 1 else n * Fact(n - 1) end end # 下降階乗冪 def FallingFact(x, n) if (n <= 1) x else x * FallingFact(x - 1, n - 1) end end # 順列 (異なる 10 個のものから 5 個取ってできる順列の総数) n = 10 r = 5 puts Fact(n) / Fact(n - r) puts FallingFact(n, r)
# 重複順列 (異なる 10 個のものから重複を許して 5 個取ってできる順列の総数) n = 10 r = 5 puts n ** r
# 組合せ def Comb(n, r) if (r == 0 || r == n) 1 elsif (r == 1) n else Comb(n - 1, r - 1) + Comb(n - 1, r) end end # 組合せ (異なる 10 個のものから 5 個取ってできる組合せの総数) n = 10 r = 5 puts Comb(n, r)
# 組合せ def Comb(n, r) if (r == 0 || r == n) 1 elsif (r == 1) n else Comb(n - 1, r - 1) + Comb(n - 1, r) end end # 重複組合せ (異なる 10 個のものから重複を許して 5 個とる組合せの総数) n = 10 r = 5 puts Comb(n + r - 1, r)
# 自作の正弦関数 def mySin(x, n, nega, numerator, denominator, y) m = 2 * n denominator = denominator * (m + 1) * m numerator = numerator * x * x a = numerator / denominator # 十分な精度になったら処理を抜ける if (a <= 0.00000000001) y else y + mySin(x, n + 1, !nega, numerator, denominator, nega ? a : -a) end end (1..24).each do |i| degree = i * 15 if degree % 30 == 0 || degree % 45 == 0 radian = degree * Math::PI / 180 # 自作の正弦関数 d1 = mySin(radian, 1, false, radian, 1.0, radian) # 標準の正弦関数 d2 = Math.sin(radian) # 標準関数との差異 printf("%3d : %13.10f - %13.10f = %13.10f\n", degree, d1, d2, d1 - d2) end end
# 自作の余弦関数 def myCos(x, n, nega, numerator, denominator, y) m = 2 * n denominator = denominator * m * (m - 1) numerator = numerator * x * x a = numerator / denominator # 十分な精度になったら処理を抜ける if (a <= 0.00000000001) y else y + myCos(x, n + 1, !nega, numerator, denominator, nega ? a : -a) end end (0..24).each do |i| degree = i * 15 if degree % 30 == 0 || degree % 45 == 0 radian = degree * Math::PI / 180 # 自作の余弦関数 d1 = myCos(radian, 1, false, 1.0, 1.0, 1.0) # 標準の余弦関数 d2 = Math.cos(radian) # 標準関数との差異 printf("%3d : %13.10f - %13.10f = %13.10f\n", degree, d1, d2, d1 - d2) end end
# 自作の正接関数 def myTan(x, x2, n, t) t = x2 / (n - t) n -= 2 if (n <= 1) x / (1 - t) else myTan(x, x2, n, t) end end (0..12).each do |i| if i * 15 % 180 != 0 degree = i * 15 - 90 radian = degree * Math::PI / 180 x2 = radian * radian # 自作の正接関数 d1 = myTan(radian, x2, 15, 0.0) # 15:必要な精度が得られる十分大きな奇数 # 標準の正接関数 d2 = Math.tan(radian) # 標準関数との差異 printf("%3d : %13.10f - %13.10f = %13.10f\n", degree, d1, d2, d1 - d2) end end
# 自作の指数関数 def myExp(x, n, numerator, denominator, y) denominator = denominator * n numerator = numerator * x a = numerator / denominator; # 十分な精度になったら処理を抜ける if (a.abs <= 0.00000000001) y else y + myExp(x, n + 1, numerator, denominator, a) end end (0..20).each do |i| x = (i - 10) / 4.0 # 標準の指数関数 d1 = Math.exp(x) # 自作の指数関数 d2 = myExp(x, 1, 1.0, 1.0, 1.0) # 標準関数との差異 printf("%5.2f : %13.10f - %13.10f = %13.10f\n", x, d1, d2, d1 - d2) end
# 自作の指数関数 def myExp(x, x2, n, t) t = x2 / (n + t) n -= 4 if (n < 6) 1 + ((2 * x) / (2 - x + t)) else myExp(x, x2, n, t) end end (0..20).each do |i| x = (i - 10) / 4.0 # 標準の指数関数 d1 = Math.exp(x) # 自作の指数関数 x2 = x * x d2 = myExp(x, x2, 30, 0.0) # 30:必要な精度が得られるよう, 6から始めて4ずつ増加させる # 標準関数との差異 printf("%5.2f : %13.10f - %13.10f = %13.10f\n", x, d1, d2, d1 - d2) end
# 自作の対数関数 def myLog(x2, numerator, denominator, y) denominator = denominator + 2 numerator = numerator * x2 * x2 a = numerator / denominator # 十分な精度になったら処理を抜ける if (a.abs <= 0.00000000001) y else y + myLog(x2, numerator, denominator, a) end end (1..20).each do |i| x = i / 5.0 # 標準の対数関数 d1 = Math.log(x) # 自作の対数関数 x2 = (x - 1) / (x + 1) d2 = 2 * myLog(x2, x2, 1.0, x2) # 標準関数との差異 printf("%5.2f : %13.10f - %13.10f = %13.10f\n", x, d1, d2, d1 - d2) end
# 自作の対数関数 def myLog(x, n, t) n2 = n x2 = x if (n > 3) if (n % 2 == 0) n2 = 2 end x2 = x * (n / 2) end t = x2 / (n2 + t) if (n <= 2) x / (1 + t); else myLog(x, n - 1, t) end end (1..20).each do |i| x = i / 5.0 # 標準の対数関数 d1 = Math.log(x) # 自作の対数関数 d2 = myLog(x - 1, 27, 0.0) # 27:必要な精度が得られる十分大きな奇数 # 標準関数との差異 printf("%5.2f : %13.10f - %13.10f = %13.10f\n", x, d1, d2, d1 - d2) end
# 自作の双曲線正弦関数 def mySinh(x, n, numerator, denominator, y) m = 2 * n denominator = denominator * (m + 1) * m numerator = numerator * x * x a = numerator / denominator # 十分な精度になったら処理を抜ける if (a.abs <= 0.00000000001) y else y + mySinh(x, n + 1, numerator, denominator, a) end end (0..20).each do |i| x = i - 10 # 自作の双曲線正弦関数 d1 = mySinh(x, 1, x, 1.0, x) # 標準の双曲線正弦関数 d2 = Math.sinh(x) # 標準関数との差異 printf("%3d : %17.10f - %17.10f = %13.10f\n", x, d1, d2, d1 - d2) end
# 自作の双曲線余弦関数 def myCosh(x, n, numerator, denominator, y) m = 2 * n denominator = denominator * m * (m - 1) numerator = numerator * x * x a = numerator / denominator # 十分な精度になったら処理を抜ける if (a.abs <= 0.00000000001) y else y + myCosh(x, n + 1, numerator, denominator, a) end end (0..20).each do |i| x = i - 10 # 自作の双曲線余弦関数 d1 = myCosh(x, 1, 1.0, 1.0, 1.0) # 標準の双曲線余弦関数 d2 = Math.cosh(x) # 標準関数との差異 printf("%3d : %17.10f - %17.10f = %13.10f\n", x, d1, d2, d1 - d2) end
def f(x) 4 / (1 + x * x) end a = 0 b = 1 # 台形則で積分 n = 2 (1..10).each do |j| h = (b - a) / n.to_f s = 0 x = a (1..(n - 1)).each do |i| x += h s += f(x) end s = h * ((f(a) + f(b)) / 2 + s) n *= 2 # 結果を π と比較 printf("%2d : %13.10f, %13.10f\n", j, s, s - Math::PI) end
def f(x) 4 / (1 + x * x) end a = 0 b = 1 # 中点則で積分 n = 2 (1..10).each do |j| h = (b - a) / n.to_f s = 0 x = a + (h / 2) (1..n).each do |i| s += f(x) x += h end s *= h n *= 2 # 結果を π と比較 printf("%2d : %13.10f, %13.10f\n", j, s, s - Math::PI) end
def f(x) 4 / (1 + x * x) end a = 0 b = 1 # Simpson則で積分 n = 2 (1..5).each do |j| h = (b - a) / n.to_f s2 = 0 s4 = 0 x = a + h (1..(n / 2)).each do |i| s4 += f(x) x += h s2 += f(x) x += h end s2 = (s2 - f(b)) * 2 + f(a) + f(b) s4 *= 4 s = (s2 + s4) * h / 3 n *= 2 # 結果を π と比較 printf("%2d : %13.10f, %13.10f\n", j, s, s - Math::PI) end
def f(x) 4 / (1 + x * x) end a = 0 b = 1 t = Array.new(7) { Array.new(7) } # 台形則で積分 n = 2 (1..6).each do |i| h = (b - a) / n.to_f s = 0 x = a (1..(n - 1)).each do |j| x += h s += f(x) end # 結果を保存 t[i][1] = h * ((f(a) + f(b)) / 2 + s) n *= 2 end # Richardsonの補外法 n = 4 (2..6).each do |j| (j..6).each do |i| t[i][j] = t[i][j - 1] + (t[i][j - 1] - t[i - 1][j - 1]) / (n - 1) if i == j # 結果を π と比較 printf("%2d : %13.10f, %13.10f\n", j, t[i][j], t[i][j] - Math::PI) end end n *= 4 end
# データ点の数 - 1 N = 6 # 元の関数 def f(x) x - (x ** 3) / (3 * 2) + (x ** 5) / (5 * 4 * 3 * 2) end # Lagrange (ラグランジュ) 補間 def lagrange(d, x, y) sum = 0 (0..N).each do |i| prod = y[i] (0..N).each do |j| if j != i prod *= (d - x[j]) / (x[i] - x[j]) end end sum += prod end sum end x = Array.new(N) y = Array.new(N) # 1.5刻みで -4.5〜4.5 まで, 7点だけ値をセット (0..N).each do |i| d = i * 1.5 - 4.5 x[i] = d y[i] = f(d) end # 0.5刻みで 与えられていない値を補間 (0..18).each do |i| d = i * 0.5 - 4.5 d1 = f(d) d2 = lagrange(d, x, y) # 元の関数と比較 printf("%5.2f\t%8.5f\t%8.5f\t%8.5f\n", d, d1, d2, d1 - d2) end
# データ点の数 - 1 N = 6 # 元の関数 def f(x) x - (x ** 3) / (3 * 2) + (x ** 5) / (5 * 4 * 3 * 2) end # Neville (ネヴィル) 補間 def neville(d, x, y) w = Array.new(N+1) { Array.new(N) } (0..N).each do |i| w[0][i] = y[i] end (1..N).each do |j| (0..(N - j)).each do |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]) end end w[N][0] end x = Array.new(N) y = Array.new(N) # 1.5刻みで -4.5〜4.5 まで, 7点だけ値をセット (0..N).each do |i| d = i * 1.5 - 4.5 x[i] = d y[i] = f(d) end # 0.5刻みで 与えられていない値を補間 (0..18).each do |i| d = i * 0.5 - 4.5 d1 = f(d) d2 = neville(d, x, y) # 元の関数と比較 printf("%5.2f\t%8.5f\t%8.5f\t%8.5f\n", d, d1, d2, d1 - d2) end
# データ点の数 - 1 N = 6 # 元の関数 def f(x) x - (x ** 3) / (3 * 2) + (x ** 5) / (5 * 4 * 3 * 2) end # Newton (ニュートン) 補間 def newton(d, x, a) sum = a[0] (1..N).each do |i| prod = a[i] (0..(i - 1)).each do |j| if j != i prod *= (d - x[j]) end end sum += prod end sum end x = Array.new(N) y = Array.new(N) # 1.5刻みで -4.5〜4.5 まで, 7点だけ値をセット (0..N).each do |i| d = i * 1.5 - 4.5 x[i] = d y[i] = f(d) end # 差分商の表を作る d = Array.new(N+1) { Array.new(N) } (0..N).each do |j| d[0][j] = y[j] end (1..N).each do |i| (0..(N-i)).each do |j| d[i][j] = (d[i-1][j+1] - d[i-1][j]) / (x[j+i] - x[j]) end end # n階差分商 a = Array.new(N) (0..N).each do |j| a[j] = d[j][0] end # 0.5刻みで 与えられていない値を補間 (0..18).each do |i| d1 = i * 0.5 - 4.5 d2 = f(d1) d3 = newton(d1, x, a) # 元の関数と比較 printf("%5.2f\t%8.5f\t%8.5f\t%8.5f\n", d1, d2, d3, d2 - d3) end
# データ点の数 - 1 N = 6 Nx2 = 13 # 元の関数 def f(x) x - (x ** 3) / (3 * 2) + (x ** 5) / (5 * 4 * 3 * 2) end # 導関数 def fd(x) 1 - (x ** 2) / 2 + (x ** 4) / (4 * 3 * 2) end # Hermite (エルミート) 補間 def hermite(d, z, a) sum = a[0] (1..Nx2).each do |i| prod = a[i] (0..(i - 1)).each do |j| if j != i prod *= (d - z[j]) end end sum += prod end sum end x = Array.new(N) y = Array.new(N) yd = Array.new(N) # 1.5刻みで -4.5〜4.5 まで, 7点だけ値をセット (0..N).each do |i| d = i * 1.5 - 4.5 x[i] = d y[i] = f(d) yd[i] = fd(d) end # 差分商の表を作る z = Array.new(Nx2) d = Array.new(Nx2+1) { Array.new(Nx2) } (0..Nx2).each do |i| j = i / 2 z[i] = x[j] d[0][i] = y[j] end (1..Nx2).each do |i| (0..(Nx2-i)).each do |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]) end end end # n階差分商 a = Array.new(Nx2) (0..Nx2).each do |j| a[j] = d[j][0] end # 0.5刻みで 与えられていない値を補間 (0..18).each do |i| d1 = i * 0.5 - 4.5 d2 = f(d1) d3 = hermite(d1, z, a) # 元の関数と比較 printf("%5.2f\t%8.5f\t%8.5f\t%8.5f\n", d1, d2, d3, d2 - d3) end
# データ点の数 - 1 N = 6 # 元の関数 def f(x) x - (x ** 3) / (3 * 2) + (x ** 5) / (5 * 4 * 3 * 2) end # Spline (スプライン) 補間 def spline(d, x, y, z) # 補間関数値がどの区間にあるか k = -1 (1..N).each do |i| if d <= x[i] k = i - 1 break end end if k < 0 k = N end d1 = x[k+1] - d d2 = d - x[k] d3 = x[k+1] - x[k] (z[k] * (d1 ** 3) + z[k+1] * (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 end x = Array.new(N) y = Array.new(N) # 1.5刻みで -4.5〜4.5 まで, 7点だけ値をセット (0..N).each do |i| d = i * 1.5 - 4.5 x[i] = d y[i] = f(d) end # 3項方程式の係数の表を作る a = Array.new(N) b = Array.new(N) c = Array.new(N) d = Array.new(N) (1..(N - 1)).each do |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])) end # 3項方程式を解く (ト−マス法) g = Array.new(N) s = Array.new(N) g[1] = b[1] s[1] = d[1] (2..(N - 1)).each do |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] end z = Array.new(N) z[0] = 0 z[N] = 0 z[N-1] = s[N-1] / g[N-1] i = N - 2 while i > 0 do z[i] = (s[i] - c[i] * z[i+1]) / g[i] i -= 1 end # 0.5刻みで 与えられていない値を補間 (0..18).each do |i| d1 = i * 0.5 - 4.5 d2 = f(d1) d3 = spline(d1, x, y, z) # 元の関数と比較 printf("%5.2f\t%8.5f\t%8.5f\t%8.5f\n", d1, d2, d3, d2 - d3) end
# 重力加速度 $g = -9.8 # 空気抵抗係数 $k = -0.01 # 時間間隔(秒) h = 0.01 # 角度 degree = 45 radian = degree * Math::PI / 180.0 # 初速 250 km/h -> 秒速に変換 v = 250 * 1000 / 3600 # 水平方向の速度 vx = Array.new(2) vx[0] = v * Math.cos(radian) # 鉛直方向の速度 vy = Array.new(2) vy[0] = v * Math.sin(radian) # 経過秒数 t = 0.0 # 位置 x = 0.0 y = 0.0 # 空気抵抗による水平方向の減速分 def fx(vx, vy) return $k * Math.sqrt(vx * vx + vy * vy) * vx end # 重力と空気抵抗による鉛直方向の減速分 def fy(vx, vy) return $g + ($k * Math.sqrt(vx * vx + vy * vy) * vy) end # Euler法 i = 1 while y >= 0.0 do # 経過秒数 t = i * h # 位置 x += h * vx[0] y += h * vy[0] 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] i += 1 end
# 重力加速度 $g = -9.8 # 空気抵抗係数 $k = -0.01 # 時間間隔(秒) h = 0.01 # 角度 degree = 45 radian = degree * Math::PI / 180.0 # 初速 250 km/h -> 秒速に変換 v = 250 * 1000 / 3600 # 水平方向の速度 vx = Array.new(3) vx[0] = v * Math.cos(radian) # 鉛直方向の速度 vy = Array.new(3) vy[0] = v * Math.sin(radian) # 経過秒数 t = 0.0 # 位置 x = Array.new(3) x[0] = 0.0 y = Array.new(3) y[0] = 0.0 # 空気抵抗による水平方向の減速分 def fx(vx, vy) return $k * Math.sqrt(vx * vx + vy * vy) * vx end # 重力と空気抵抗による鉛直方向の減速分 def fy(vx, vy) return $g + ($k * Math.sqrt(vx * vx + vy * vy) * vy) end # Heun法 i = 1 while y[0] >= 0.0 do # 経過秒数 t = 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] printf("%4.2f\t%8.5f\t%9.5f\t%9.5f\t%8.5f\n", t, vx[0], vy[0], x[0], y[0]) i += 1 end
# 重力加速度 $g = -9.8 # 空気抵抗係数 $k = -0.01 # 時間間隔(秒) h = 0.01 # 角度 degree = 45 radian = degree * Math::PI / 180.0 # 初速 250 km/h -> 秒速に変換 v = 250 * 1000 / 3600 # 水平方向の速度 vx = Array.new(2) vx[0] = v * Math.cos(radian) # 鉛直方向の速度 vy = Array.new(2) vy[0] = v * Math.sin(radian) # 経過秒数 t = 0.0 # 位置 x = Array.new(2) x[0] = 0.0 y = Array.new(2) y[0] = 0.0 # 空気抵抗による水平方向の減速分 def fx(vx, vy) return $k * Math.sqrt(vx * vx + vy * vy) * vx end # 重力と空気抵抗による鉛直方向の減速分 def fy(vx, vy) return $g + ($k * Math.sqrt(vx * vx + vy * vy) * vy) end # 中点法 i = 1 while y[0] >= 0.0 do # 経過秒数 t = i * h # 位置・速度 vx[1] = h * fx(vx[0], vy[0]) vy[1] = h * fy(vx[0], vy[0]) wx = vx[0] + vx[1] / 2.0 wy = vy[0] + vy[1] / 2.0 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 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]) i += 1 end
# 重力加速度 $g = -9.8 # 空気抵抗係数 $k = -0.01 # 時間間隔(秒) h = 0.01 # 角度 degree = 45 radian = degree * Math::PI / 180.0 # 初速 250 km/h -> 秒速に変換 v = 250 * 1000 / 3600 # 水平方向の速度 vx = Array.new(5) vx[0] = v * Math.cos(radian) # 鉛直方向の速度 vy = Array.new(5) vy[0] = v * Math.sin(radian) # 経過秒数 t = 0.0 # 位置 x = Array.new(5) x[0] = 0.0 y = Array.new(5) y[0] = 0.0 # 空気抵抗による水平方向の減速分 def fx(vx, vy) return $k * Math.sqrt(vx * vx + vy * vy) * vx end # 重力と空気抵抗による鉛直方向の減速分 def fy(vx, vy) return $g + ($k * Math.sqrt(vx * vx + vy * vy) * vy) end # Runge-Kutta法 i = 1 while y[0] >= 0.0 do # 経過秒数 t = 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]) wx = vx[0] + vx[1] / 2 wy = 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 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]) i += 1 end
# 重力加速度 $g = -9.8 # 空気抵抗係数 $k = -0.01 # 時間間隔(秒) h = 0.01 # 角度 degree = 45 radian = degree * Math::PI / 180.0 # 初速 250 km/h -> 秒速に変換 v = 250 * 1000 / 3600 # 水平方向の速度 vx = Array.new(5) vx[0] = v * Math.cos(radian) # 鉛直方向の速度 vy = Array.new(5) vy[0] = v * Math.sin(radian) # 経過秒数 t = 0.0 # 位置 x = Array.new(5) x[0] = 0.0 y = Array.new(5) y[0] = 0.0 # 空気抵抗による水平方向の減速分 def fx(vx, vy) return $k * Math.sqrt(vx * vx + vy * vy) * vx end # 重力と空気抵抗による鉛直方向の減速分 def fy(vx, vy) return $g + ($k * Math.sqrt(vx * vx + vy * vy) * vy) end # Runge-Kutta-Gill法 i = 1 while y[0] >= 0.0 do # 経過秒数 t = 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]) wx = vx[0] + vx[1] / 2 wy = 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 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]) i += 1 end
include Math def f(x) x * x - 2.0 end 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 - sqrt(2)) fc = f(c) if fc.abs < 0.0000000001 break end if fc < 0 # f(c) < 0 であれば, 解は区間 (c, b) の中に存在 a = c else # f(c) > 0 であれば, 解は区間 (a, c) の中に存在 b = c end end c end a = 1.0 b = 2.0 printf("%12.10f\n", bisection(a, b))
include Math def f(x) x * x - 2.0 end 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 - sqrt(2)) fc = f(c) if fc.abs < 0.0000000001 break end if fc < 0 # f(c) < 0 であれば, 解は区間 (c, b) の中に存在 a = c else # f(c) > 0 であれば, 解は区間 (a, c) の中に存在 b = c end end c end a = 1.0 b = 2.0 printf("%12.10f\n", falseposition(a, b))
include Math def g(x) (x / 2.0) + (1.0 / x) end def iterative(x0) while true x1 = g(x0) printf("%12.10f\t%13.10f\n", x1, x1 - sqrt(2)) if (x1 - x0).abs < 0.0000000001 break end x0 = x1 end x1 end x = 1.0 printf("%12.10f\n", iterative(x))
include Math def f0(x) x * x - 2.0 end def f1(x) 2.0 * x end def newton(x0) while true x1 = x0 - (f0(x0) / f1(x0)) printf("%12.10f\t%13.10f\n", x1, x1 - sqrt(2)) if (x1 - x0).abs < 0.0000000001 break end x0 = x1 end x1 end x = 2.0 printf("%12.10f\n", newton(x))
include Math def f0(x) x * x - 2.0 end def f1(x) 2.0 * x end def f2(x) 2.0 end 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 - sqrt(2)) if (x1 - x0).abs < 0.0000000001 break end x0 = x1 end x1 end x = 2.0 printf("%12.10f\n", bailey(x))
include Math def f(x) x * x - 2.0 end 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 - sqrt(2)) if (x2 - x1).abs < 0.0000000001 break end x0 = x1 x1 = x2 end x2 end x0 = 1.0 x1 = 2.0 printf("%12.10f\n", secant(x0, x1))
# coding: Shift_JIS N = 3 # 1次元配列を表示 def disp_vector(row) row.each do |col| printf("%14.10f\t", col) end puts "" end # ヤコビの反復法 def jacobi(a, b, x0) while true x1 = [] finish = true (0..N).each do |i| x1[i] = 0 (0..N).each do |j| if (j != i) x1[i] += a[i][j] * x0[j] end end x1[i] = (b[i] - x1[i]) / a[i][i] if ((x1[i] - x0[i]).abs > 0.0000000001) finish = false end end (0..N).each do |i| x0[i] = x1[i] end return if finish disp_vector(x0) end end 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]] b = [20.0, 16.0, 8.0, 17.0] c = [ 0.0, 0.0, 0.0, 0.0] # ヤコビの反復法 jacobi(a,b,c) puts "X" disp_vector(c)
# coding: Shift_JIS N = 3 # 1次元配列を表示 def disp_vector(row) row.each do |col| printf("%14.10f\t", col) end puts "" end # ガウス・ザイデル法 def gauss(a, b, x0) while true finish = true (0..N).each do |i| x1 = 0 (0..N).each do |j| if (j != i) x1 += a[i][j] * x0[j] end end x1 = (b[i] - x1) / a[i][i] if ((x1 - x0[i]).abs > 0.0000000001) finish = false end x0[i] = x1 end return if (finish) disp_vector(x0) end end 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]] b = [20.0, 16.0, 8.0, 17.0] c = [ 0.0, 0.0, 0.0, 0.0] # ガウス・ザイデル法 gauss(a,b,c) puts "X" disp_vector(c)
# coding: Shift_JIS N = 3 # 1次元配列を表示 def disp_matrix(matrix) matrix.each do |row| row.each do |col| printf("%14.10f\t", col) end puts "" end end # 2次元配列を表示 def disp_vector(row) row.each do |col| printf("%14.10f\t", col) end puts "" end # 前進消去 def forward_elimination(a, b) 0.upto(N - 1) do |pivot| (pivot + 1).upto(N) do |row| s = a[row][pivot] / a[pivot][pivot] pivot.upto(N) do |col| a[row][col] -= a[pivot][col] * s end b[row] -= b[pivot] * s end end end # 後退代入 def backward_substitution(a, b) N.downto(0) do |row| N.downto(row + 1) do |col| b[row] -= a[row][col] * b[col] end b[row] /= a[row][row] end end # ピボット選択 def pivoting(a, b) (0..N).each do |pivot| # 各列で 一番値が大きい行を 探す max_row = pivot max_val = 0 (pivot..N).each do |row| if ((a[row][pivot]).abs > max_val) # 一番値が大きい行 max_val = (a[row][pivot]).abs max_row = row end end # 一番値が大きい行と入れ替え if (max_row != pivot) tmp = 0 (0..N).each do |col| tmp = a[max_row][col] a[max_row][col] = a[pivot][col] a[pivot][col] = tmp end tmp = b[max_row] b[max_row] = b[pivot] b[pivot] = tmp end end end 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]] b = [8.0,17.0,20.0,16.0] # ピボット選択 pivoting(a,b) puts "pivoting" puts "A" disp_matrix(a) puts "B" disp_vector(b) puts "" # 前進消去 forward_elimination(a,b) puts "forward elimination" puts "A" disp_matrix(a) puts "B" disp_vector(b) puts "" # 後退代入 backward_substitution(a,b) puts "X" disp_vector(b)
# coding: Shift_JIS N = 3 # 2次元配列を表示 def disp_matrix(matrix) matrix.each do |row| row.each do |col| printf("%14.10f\t", col) end puts "" end end # 1次元配列を表示 def disp_vector(row) row.each do |col| printf("%14.10f\t", col) end puts "" end # 前進消去 def forward_elimination(a, b) 0.upto(N) do |pivot| 0.upto(N) do |row| next if (row == pivot) s = a[row][pivot] / a[pivot][pivot] pivot.upto(N) do |col| a[row][col] -= a[pivot][col] * s end b[row] -= b[pivot] * s end end end # 後退代入 def backward_substitution(a, b) 0.upto(N) do |pivot| b[pivot] /= a[pivot][pivot] end end # ピボット選択 def pivoting(a, b) (0..N).each do |pivot| # 各列で 一番値が大きい行を 探す max_row = pivot max_val = 0 (pivot..N).each do |row| if ((a[row][pivot]).abs > max_val) # 一番値が大きい行 max_val = (a[row][pivot]).abs max_row = row end end # 一番値が大きい行と入れ替え if (max_row != pivot) tmp = 0 (0..N).each do |col| tmp = a[max_row][col] a[max_row][col] = a[pivot][col] a[pivot][col] = tmp end tmp = b[max_row] b[max_row] = b[pivot] b[pivot] = tmp end end end # 状態を確認 def disp_progress(a, b) (0..N).each do |i| (0..N).each do |j| printf("%14.10f\t", a[i][j]) end printf("%14.10f\n", b[i]) end puts "" end 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]] b = [8.0,17.0,20.0,16.0] # ピボット選択 pivoting(a,b) puts "pivoting" puts "A" disp_matrix(a) puts "B" disp_vector(b) puts "" # 前進消去 forward_elimination(a,b) puts "forward elimination" puts "A" disp_matrix(a) puts "B" disp_vector(b) puts "" # 後退代入 backward_substitution(a,b) puts "X" disp_vector(b)
# coding: Shift_JIS N = 3 # 2次元配列を表示 def disp_matrix(matrix) matrix.each do |row| row.each do |col| printf("%14.10f\t", col) end puts "" end end # 1次元配列を表示 def disp_vector(row) row.each do |col| printf("%14.10f\t", col) end puts "" end # 前進消去 def forward_elimination(a, b) 0.upto(N - 1) do |pivot| (pivot + 1).upto(N) do |row| s = a[row][pivot] / a[pivot][pivot] pivot.upto(N) do |col| a[row][col] -= a[pivot][col] * s # これが 上三角行列 end a[row][pivot] = s # これが 下三角行列 # b[row] -= b[pivot] * s # この値は変更しない end end end # Ly=b から y を求める (前進代入) def forward_substitution(a, b, y) (0..N).each do |row| (0..row).each do |col| b[row] -= a[row][col] * y[col] end y[row] = b[row] end end # Ux=y から x を求める (後退代入) def backward_substitution(a, y, x) N.downto(0) do |row| N.downto(row + 1) do |col| y[row] -= a[row][col] * x[col] end x[row] = y[row] / a[row][row] end end # ピボット選択 def pivoting(a, b) (0..N).each do |pivot| # 各列で 一番値が大きい行を 探す max_row = pivot max_val = 0 (pivot..N).each do |row| if ((a[row][pivot]).abs > max_val) # 一番値が大きい行 max_val = (a[row][pivot]).abs max_row = row end end # 一番値が大きい行と入れ替え if (max_row != pivot) tmp = 0 (0..N).each do |col| tmp = a[max_row][col] a[max_row][col] = a[pivot][col] a[pivot][col] = tmp end tmp = b[max_row] b[max_row] = b[pivot] b[pivot] = tmp end end end 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]] b = [8.0,17.0,20.0,16.0] # ピボット選択 pivoting(a,b) puts "pivoting" puts "A" disp_matrix(a) puts "B" disp_vector(b) puts "" # 前進消去 forward_elimination(a,b) puts "LU" disp_matrix(a) # Ly=b から y を求める (前進代入) y = [0.0,0.0,0.0,0.0] forward_substitution(a,b,y) puts "Y" disp_vector(y) # Ux=y から x を求める (後退代入) x = [0.0,0.0,0.0,0.0] backward_substitution(a,y,x) puts "X" disp_vector(x)
# coding: Shift_JIS N = 3 # 2次元配列を表示 def disp_matrix(matrix) matrix.each do |row| row.each do |col| printf("%14.10f\t", col) end puts "" end end # 1次元配列を表示 def disp_vector(row) row.each do |col| printf("%14.10f\t", col) end puts "" end # 前進消去 def forward_elimination(a, b) 0.upto(N) do |pivot| s = 0 0.upto(pivot-1) do |col| s += a[pivot][col] * a[pivot][col] end # ここで根号の中が負の値になると計算できない! a[pivot][pivot] = Math.sqrt(a[pivot][pivot] - s) (pivot + 1).upto(N) do |row| s = 0 0.upto(pivot-1) do |col| s += a[row][col] * a[pivot][col] end a[row][pivot] = (a[row][pivot] - s) / a[pivot][pivot] a[pivot][row] = a[row][pivot] end end end # 前進代入 def forward_substitution(a, b, y) (0..N).each do |row| (0..row).each do |col| b[row] -= a[row][col] * y[col] end y[row] = b[row] / a[row][row] end end # 後退代入 def backward_substitution(a, y, x) N.downto(0) do |row| N.downto(row + 1) do |col| y[row] -= a[row][col] * x[col] end x[row] = y[row] / a[row][row] end end 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]] b = [34.0,68.0,96.0,125.0] puts "A" disp_matrix(a) puts "B" disp_vector(b) puts "" # 前進消去 forward_elimination(a,b) puts "LL^T" disp_matrix(a) # Ly=b から y を求める (前進代入) y = [0.0,0.0,0.0,0.0] forward_substitution(a,b,y) puts "Y" disp_vector(y) # L^Tx=y から x を求める (後退代入) x = [0.0,0.0,0.0,0.0] backward_substitution(a,y,x) puts "X" disp_vector(x)
# coding: Shift_JIS N = 3 # 2次元配列を表示 def disp_matrix(matrix) matrix.each do |row| row.each do |col| printf("%14.10f\t", col) end puts "" end end # 1次元配列を表示 def disp_vector(row) row.each do |col| printf("%14.10f\t", col) end puts "" end # 前進消去 def forward_elimination(a, b) (0..N).each do |pivot| # pivot < k の場合 s = 0 (0..(pivot-1)).each do |col| s = a[pivot][col] (0..(col-1)).each do |k| s -= a[pivot][k] * a[col][k] * a[k][k] end a[pivot][col] = s / a[col][col] a[col][pivot] = a[pivot][col] end # pivot == k の場合 s = a[pivot][pivot] (0..(pivot-1)).each do |k| s -= a[pivot][k] * a[pivot][k] * a[k][k] end a[pivot][pivot] = s end end # 前進代入 def forward_substitution(a, b, y) (0..N).each do |row| (0..row).each do |col| b[row] -= a[row][col] * y[col] end y[row] = b[row] end end # 後退代入 def backward_substitution(a, y, x) N.downto(0) do |row| N.downto(row + 1) do |col| y[row] -= a[row][col] * a[row][row] * x[col] end x[row] = y[row] / a[row][row] end end 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]] b = [34.0,68.0,96.0,125.0] puts "A" disp_matrix(a) puts "B" disp_vector(b) puts "" # 前進消去 forward_elimination(a,b) puts "LDL^T" disp_matrix(a) # Ly=b から y を求める (前進代入) y = [0.0,0.0,0.0,0.0] forward_substitution(a,b,y) puts "Y" disp_vector(y) # DL^Tx=y から x を求める (後退代入) x = [0.0,0.0,0.0,0.0] backward_substitution(a,y,x) puts "X" disp_vector(x)
# coding: Shift_JIS N = 3 # 1次元配列を表示 def disp_vector(row) row.each do |col| printf("%14.10f\t", col) end puts "" end # 正規化 (ベクトルの長さを1にする) def normarize(x) s = 0.0 (0..N).each do |i| s += x[i] * x[i] end s = Math.sqrt(s) (0..N).each do |i| x[i] /= s end end # ベキ乗法 def power(a, x0) lambda = 0.0 # 正規化 (ベクトル x0 の長さを1にする) normarize(x0) e0 = 0.0 (0..N).each do |i| e0 += x0[i] end (1..200).each do |k| # 1次元配列を表示 printf("%3d\t", k) disp_vector(x0) # 行列の積 x1 = A × x0 x1 = [ 0.0, 0.0, 0.0, 0.0] (0..N).each do |i| (0..N).each do |j| x1[i] += a[i][j] * x0[j] end end # 内積 p0 = 0.0 p1 = 0.0 (0..N).each do |i| p0 += x1[i] * x1[i] p1 += x1[i] * x0[i] end # 固有値 lambda = p0 / p1 # 正規化 (ベクトル x1 の長さを1にする) normarize(x1) # 収束判定 e1 = 0.0 (0..N).each do |i| e1 += x1[i] end break if ((e0 - e1).abs < 0.00000000001) (0..N).each do |i| x0[i] = x1[i] end e0 = e1 end lambda end # ベキ乗法で最大固有値を求める 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]] x = [1.0 ,0.0 ,0.0 ,0.0] # ベキ乗法 lambda0 = power(a, x) puts "" puts "eigenvalue" printf("%14.10f\n", lambda0) puts "eigenvector" disp_vector(x)
# coding: Shift_JIS N = 3 # 前進消去 def forward_elimination(a) 0.upto(N - 1) do |pivot| (pivot + 1).upto(N) do |row| s = a[row][pivot] / a[pivot][pivot] pivot.upto(N) do |col| a[row][col] -= a[pivot][col] * s # これが 上三角行列 end a[row][pivot] = s # これが 下三角行列 end end end # 前進代入 def forward_substitution(a, y, b) (0..N).each do |row| (0..row).each do |col| b[row] -= a[row][col] * y[col] end y[row] = b[row] end end # 後退代入 def backward_substitution(a, x, y) (0..N).reverse_each do |row| ((row + 1)..N).reverse_each do |col| y[row] -= a[row][col] * x[col] end x[row] = y[row] / a[row][row] end end # 1次元配列を表示 def disp_vector(row) row.each do |col| printf("%14.10f\t", col) end puts "" end # 正規化 (ベクトルの長さを1にする) def normarize(x) s = 0.0 (0..N).each do |i| s += x[i] * x[i] end s = Math.sqrt(s) (0..N).each do |i| x[i] /= s end end # 逆ベキ乗法 def inverse(a, x0) lambda = 0.0 # 正規化 (ベクトル x0 の長さを1にする) normarize(x0) e0 = 0.0 (0..N).each do |i| e0 += x0[i] end (1..200).each do |k| # 1次元配列を表示 printf("%3d\t", k) disp_vector(x0) # Ly=b から y を求める (前進代入) y = [0.0, 0.0, 0.0, 0.0] b = [0.0, 0.0, 0.0, 0.0] (0..N).each do |i| b[i] = x0[i] end 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 (0..N).each do |i| p0 += x1[i] * x1[i] p1 += x1[i] * x0[i] end # 固有値 lambda = p1 / p0 # 正規化 (ベクトル x1 の長さを1にする) normarize(x1) # 収束判定 e1 = 0.0 (0..N).each do |i| e1 += x1[i] end break if ((e0 - e1).abs < 0.00000000001) (0..N).each do |i| x0[i] = x1[i] end e0 = e1 end lambda end # 逆ベキ乗法で最小固有値を求める 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]] x = [1.0 ,0.0 ,0.0 ,0.0] # LU分解 forward_elimination(a) # 逆ベキ乗法 lambda0 = inverse(a, x) puts "" puts "eigenvalue" printf("%14.10f\n", lambda0) puts "eigenvector" disp_vector(x)
# coding: Shift_JIS N = 3 # 対角要素を表示 def disp_eigenvalue(a) (0..N).each do |i| printf("%14.10f\t", a[i][i]) end puts "" end # 行列の積 def multiply(a, b, c) (0..N).each do |i| (0..N).each do |j| s = 0.0 (0..N).each do |k| s += a[i][k] * b[k][j] end c[i][j] = s end end end # LU分解 def decomp(a, l, u) (0..N).each do |i| (0..N).each do |j| l[i][j] = 0.0 u[i][j] = 0.0 end end l[0][0] = 1.0 (0..N).each do |j| u[0][j] = a[0][j] end (1..N).each do |i| u[i][0] = 0.0 l[0][i] = 0.0 l[i][0] = a[i][0] / u[0][0] end (1..N).each do |i| l[i][i] = 1.0 t = a[i][i] (0..i).each do |k| t -= l[i][k] * u[k][i] end u[i][i] = t ((i + 1)..N).each do |j| u[j][i] = 0.0 l[i][j] = 0.0 t = a[j][i] (0..i).each do |k| t -= l[j][k] * u[k][i] end l[j][i] = t / u[i][i] t = a[i][j] (0..i).each do |k| t -= l[i][k] * u[k][j] end u[i][j] = t end end end # LR分解で固有値を求める 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]] 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]] 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]] (1..200).each do |k| # LU分解 decomp(a, l, u) # 行列の積 multiply(u, l, a) # 対角要素を表示 printf("%3d\t", k) disp_eigenvalue(a) # 収束判定 e = 0.0 (1..N).each do |i| (0..(i - 1)).each do |j| e += (a[i][j]).abs end end break if (e < 0.00000000001) end puts "" puts "eigenvalue" disp_eigenvalue(a)
# coding: Shift_JIS N = 3 # 対角要素を表示 def disp_eigenvalue(a) (0..N).each do |i| printf("%14.10f\t", a[i][i]) end puts "" end # 行列の積 def multiply(a, b, c) (0..N).each do |i| (0..N).each do |j| s = 0.0 (0..N).each do |k| s += a[i][k] * b[k][j] end c[i][j] = s end end end # QR分解 def decomp(a, q, r) x = [0.0 ,0.0 ,0.0 ,0.0] (0..N).each do |k| (0..N).each do |i| x[i] = a[i][k] end (0..(k - 1)).each do |j| t = 0.0 (0..N).each do |i| t += a[i][k] * q[i][j] end r[j][k] = t r[k][j] = 0.0 (0..N).each do |i| x[i] -= t * q[i][j] end end s = 0.0 (0..N).each do |i| s += x[i] * x[i] end r[k][k] = Math.sqrt(s) (0..N).each do |i| q[i][k] = x[i] / r[k][k] end end end # QR分解で固有値を求める 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]] 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]] 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]] (1..200).each do |k| # QR分解 decomp(a, q, r) # 行列の積 multiply(r, q, a) # 対角要素を表示 printf("%3d\t", k) disp_eigenvalue(a) # 収束判定 e = 0.0 (1..N).each do |i| (0..(i - 1)).each do |j| e += (a[i][j]).abs end end break if (e < 0.00000000001) end puts "" puts "eigenvalue" disp_eigenvalue(a)
# coding: Shift_JIS N = 3 # 対角要素を表示 def disp_eigenvalue(a) (0..N).each do |i| printf("%14.10f\t", a[i][i]) end puts "" end # 1次元配列を表示 def disp_vector(row) row.each do |col| printf("%14.10f\t", col) end puts "" end # 正規化 (ベクトルの長さを1にする) def normarize(x) s = 0.0 (0..N).each do |i| s += x[i] * x[i] end s = Math.sqrt(s) (0..N).each do |i| x[i] /= s end end # 固有ベクトルを表示 def disp_eigenvector(matrix) (0..N).each do |i| row = [0.0 ,0.0 ,0.0 ,0.0] (0..N).each do |j| row[j] = matrix[i][j] end normarize(row) disp_vector(row) end end # ヤコビ法 def jacobi(a, v) (1..200).each do |k| # 最大値を探す p = 0 q = 0 max_val = 0.0 (0..N).each do |i| ((i + 1)..N).each do |j| if (max_val < (a[i][j]).abs) max_val = (a[i][j]).abs p = i q = j end end end # θ を求める t = 0.0 if ((a[p][p] - a[q][q]).abs < 0.00000000001) # a_{pp} = a_{qq} のとき、回転角tをπ/4にする t = Math::PI / 4.0 t = -t if (a[p][p] < 0) else # a_{pp} ≠ a_{qq} のとき t = Math.atan(2.0 * a[p][q] / (a[p][p] - a[q][q])) / 2.0 end # θ を使って 行列 U を作成し、A = U^t × A × U c = Math.cos(t) s = Math.sin(t) # U^t × A t1 = 0.0 t2 = 0.0 (0..N).each do |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 end # A × U (0..N).each do |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 end # 対角要素を表示 printf("%3d\t", k) disp_eigenvalue(a) # 収束判定 break if (max_val < 0.00000000001) end end # ヤコビ法で固有値を求める 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]] 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]] # ヤコビ法 jacobi(a, v) puts "" puts "eigenvalue" disp_eigenvalue(a) puts "" puts "eigenvector" disp_eigenvector(v)
# coding: Shift_JIS N = 3 # 1次元配列を表示 def disp_vector(row) row.each do |col| printf("%14.10f\t", col) end puts "" end # 2次元配列を表示 def disp_matrix(matrix) matrix.each do |row| row.each do |col| printf("%14.10f\t", col) end puts "" end end # ハウスホルダー変換 def tridiagonalize(a, d, e) v = [0.0 ,0.0 ,0.0 ,0.0] (0..(N - 2)).each do |k| w = [0.0 ,0.0 ,0.0 ,0.0] d[k] = a[k][k] t = 0.0 ((k + 1)..N).each do |i| w[i] = a[i][k] t += w[i] * w[i] end s = Math.sqrt(t) s = -s if (w[k + 1] < 0) if (s.abs < 0.00000000001) e[k + 1] = 0.0 else e[k + 1] = -s w[k + 1] += s s = 1 / Math.sqrt(w[k + 1] * s) ((k + 1)..N).each do |i| w[i] *= s end ((k + 1)..N).each do |i| s = 0.0 ((k + 1)..N).each do |j| if (j <= i) s += a[i][j] * w[j] else s += a[j][i] * w[j] end end v[i] = s end s = 0.0 ((k + 1)..N).each do |i| s += w[i] * v[i] end s /= 2.0 ((k + 1)..N).each do |i| v[i] -= s * w[i] end ((k + 1)..N).each do |i| ((k + 1)..i).each do |j| a[i][j] -= w[i] * v[j] + w[j] * v[i] end end # 次の行は固有ベクトルを求めないなら不要 ((k + 1)..N).each do |i| a[i][k] = w[i] end end end d[N - 1] = a[N - 1][N - 1] d[N] = a[N][N] e[0] = 0.0 e[N] = a[N][N - 1] # 次の行は固有ベクトルを求めないなら不要 (0..N).reverse_each do |k| w = [0.0 ,0.0 ,0.0 ,0.0] if (k < N - 1) ((k + 1)..N).each do |i| w[i] = a[i][k] end ((k + 1)..N).each do |i| s = 0.0 ((k + 1)..N).each do |j| s += a[i][j] * w[j] end v[i] = s end ((k + 1)..N).each do |i| ((k + 1)..N).each do |j| a[i][j] -= v[i] * w[j] end end end (0..N).each do |i| a[i][k] = 0.0 end a[k][k] = 1.0 end end # QR分解 def decomp(a, d, e) e[0] = 1.0 h = N - 1 h -= 1 while ((e[h]).abs < 0.00000000001) while (h > 0) e[0] = 0.0 l = h - 1 l -= 1 while ((e[l]).abs >= 0.00000000001) (1..100).each do |j| w = (d[h - 1] - d[h]) / 2.0 s = Math.sqrt(w * w + e[h] * e[h]) s = -s if (w < 0.0) x = d[l] - d[h] + e[h] * e[h] / (w + s) y = e[l + 1] z = 0.0 t = 0.0 u = 0.0 (l..(h - 1)).each do |k| if (x.abs >= y.abs) 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) s = -s if (t < 0) u = t * s end 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 # 次の行は固有ベクトルを求めないなら不要 (0..N).each do |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 end if (k < N - 1) x = e[k + 1] y = -s * e[k + 2] z = y e[k + 2] = u * e[k + 2] end end printf("%3d\t", j) disp_vector(d) # 収束判定 break if ((e[h]).abs < 0.00000000001) end e[0] = 1.0 h -= 1 while ((e[h]).abs < 0.00000000001) end # 次の行は固有ベクトルを求めないなら不要 (0..(N - 1)).each do |k| l = k ((k + 1)..N).each do |i| l = i if (d[i] > d[l]) end t = d[k] d[k] = d[l] d[l] = t (0..N).each do |i| t = a[k][i] a[k][i] = a[l][i] a[l][i] = t end end end # ハウスホルダー変換とQR分解で固有値を求める 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]] d = [1.0 ,0.0 ,0.0 ,0.0] e = [1.0 ,0.0 ,0.0 ,0.0] # ハウスホルダー変換 tridiagonalize(a, d, e) # QR分解 decomp(a, d, e) puts "" puts "eigenvalue" disp_vector(d) puts "" puts "eigenvector" disp_matrix(a)