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 までの 自然数の２乗の和
n = 1000
print (n * (n + 1) * (2 * n + 1) / 6)
# 100 までの 自然数の３乗の和
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 まで, ７点だけ値をセット
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 まで, ７点だけ値をセット
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 まで, ７点だけ値をセット
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])
# ｎ階差分商
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 まで, ７点だけ値をセット
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])
# ｎ階差分商
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 まで, ７点だけ値をセット
for i in range(0, N):
d = i * 1.5 - 4.5
x[i] = d
y[i] = f(d)
# ３項方程式の係数の表を作る
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]))
# ３項方程式を解く (ト−マス法)
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
# ヤコビの反復法
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]
print "%14.10f\t" % x0[i],
print ""
if (finish):
return
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 "解"
for i in range(0, N, 1):
print "%14.10f\t" % c[i],
print ""
# coding: Shift_JIS
N = 4
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
for i in range(0, N, 1):
print "%14.10f\t" % x0[i],
print ""
if (finish):
return
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"
for i in range(0, N, 1):
print "%14.10f\t" % c[i],
print ""
# 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
def disp_progress(a, b):
for i in range(0, N, 1):
for j in range(0, N, 1):
print "%14.10f\t" % a[i][j],
print "%14.10f" % b[i]
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"
disp_progress(a,b)
forward_elimination(a,b)
print "forward_elimination"
disp_progress(a,b)
backward_substitution(a,b)
for i in range(0, N, 1):
print "%14.10f\t" % b[i],
print ""
# coding: Shift_JIS
N = 4
def disp_progress(a, b):
for i in range(0, N, 1):
for j in range(0, N, 1):
print "%14.10f\t" % a[i][j],
print "%14.10f" % b[i]
print ""
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
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"
disp_progress(a,b)
forward_elimination(a,b)
print "forward_elimination"
disp_progress(a,b)
backward_substitution(a,b)
for i in range(0, N, 1):
print "%14.10f\t" % b[i],
print ""
# coding: Shift_JIS
N = 4
# LU
def decomp(a, b, x):
#
forward_elimination(a,b)
#
print "L"
for i in range(0, N, 1):
for j in range (0, N, 1):
if (i > j):
print "%14.10f\t" % a[i][j],
elif (i == j):
print "%14.10f\t" % 1.0,
else:
print "%14.10f\t" % 0.0,
print ""
print ""
#
print "U"
for i in range(0, N, 1):
for j in range(0, N, 1):
if (i <= j):
print "%14.10f\t" % a[i][j],
else:
print "%14.10f\t" % 0.0,
print ""
print ""
# Ly=b
y = [0.0, 0.0, 0.0, 0.0]
forward_substitution(a,b,y)
# y
print "Y"
for i in range(0, N, 1):
print "%14.10f\t" % y[i],
print ""
# Ux=y,
backward_substitution(a,y,x)
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
# Ly=b
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
def disp_progress(a, b):
for i in range(0, N, 1):
for j in range(0, N, 1):
print "%14.10f\t" % a[i][j],
print "%14.10f" % b[i]
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"
disp_progress(a,b)
x = [ 0.0, 0.0, 0.0, 0.0]
decomp(a,b,x)
print "X"
for i in range(0, N, 1):
print "%14.10f\t" % x[i],
print ""
# coding: Shift_JIS
import math
N = 4
# LL^T
def decomp(a, b, x):
#
forward_elimination(a,b)
#
print "L L^T"
for i in range(0, N, 1):
for j in range (0, N, 1):
if (j <= i):
print "%14.10f\t" % a[i][j],
else:
print "%14.10f\t" % a[j][i],
print ""
print ""
# Ly=b
y = [0.0, 0.0, 0.0, 0.0]
forward_substitution(a,b,y)
# y
print "Y"
for i in range(0, N, 1):
print "%14.10f\t" % y[i],
print ""
# Ux=y,
backward_substitution(a,y,x)
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]
# Ly=b
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[col][row] * x[col]
x[row] = y[row] / a[row][row]
def disp_progress(a):
for i in range(0, N, 1):
for j in range(0, N, 1):
print "%14.10f\t" % a[i][j],
print ""
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_progress(a)
x = [ 0.0, 0.0, 0.0, 0.0]
decomp(a,b,x)
print "X"
for i in range(0, N, 1):
print "%14.10f\t" % x[i],
print ""
# coding: Shift_JIS
N = 4
# LDL^T
def decomp(a, b, x):
#
forward_elimination(a,b)
#
print "L D"
for i in range(0, N, 1):
for j in range (0, N, 1):
if (j <= i):
print "%14.10f\t" % a[i][j],
else:
print "%14.10f\t" % a[j][i],
print ""
print ""
# Ly=b
y = [0.0, 0.0, 0.0, 0.0]
forward_substitution(a,b,y)
# y
print "Y"
for i in range(0, N, 1):
print "%14.10f\t" % y[i],
print ""
# Ux=y,
backward_substitution(a,y,x)
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]
# 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
# Ly=b
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]
def disp_progress(a):
for i in range(0, N, 1):
for j in range(0, N, 1):
print "%14.10f\t" % a[i][j],
print ""
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_progress(a)
x = [ 0.0, 0.0, 0.0, 0.0]
decomp(a,b,x)
print "X"
for i in range(0, N, 1):
print "%14.10f\t" % x[i],
print ""
|