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