Re: [資訊] python科學計算
※ 引述《gaviniscool (spectator)》之銘言：
: 最近打算用PYTHON玩玩看科學模擬
: 發現對岸好像有出電子書可以免費下載
: "用PYTHON做科學計算"
: http://hyry.dip.jp:8000/pydoc/index.html
: 最近好像有改成實體書出版
: 給有需要的朋友參考
參考書中的說明
以下是重力場三體問題的一個簡單模擬 :)
# -*- coding: utf-8 -*-
import numpy as np
import scipy as sp
import matplotlib as mpl
import matplotlib.pyplot as plt
# 為了求ODE數值解
from scipy.integrate import odeint
G = 4
pi = sp.pi
def vector(P1, P2):
x1, y1 = P1
x2, y2 = P2
return x2-x1, y2-y1
def norm(dx, dy):
return (dx**2+dy**2)**0.5
def GravityField(m1, P1, P2):
vx, vy = vector(P1, P2)
r = norm(vx, vy)
return -1*G*m1*vx*(1.0)/r**3, -1*G*m1*vy*(1.0)/r**3
def polarCoordinate(r, theta):
return r*sp.cos(2*pi*theta/360.0), r*sp.sin(2*pi*theta/360.0)
def circleOrbitSpeed(M, r):
return (G*M*1.0/r)**0.5
def tangentVector(v, theta):
return -1*v*sp.sin(2*pi*theta/360.0), v*sp.cos(2*pi*theta/360.0)
def K_Energy(m, vx, vy):
return 0.5*m*(vx**2+vy**2)
def U_Energy(m1, m2, r):
return -1*G*m1*m2/r
def ode(P, t, Ps, ms, m1, m2):
# given position vector, and parameter
# calculate dx/dt, dy/dt, dz/dt
x1, y1, vx1, vy1, x2, y2, vx2, vy2 = P # 太陽是固定的
P1 = (x1, y1)
P2 = (x2, y2)
gs1x, gs1y = GravityField(ms, Ps, P1)
gs2x, gs2y = GravityField(ms, Ps, P2)
g12x, g12y = GravityField(m1, P1, P2)
g21x, g21y = GravityField(m2, P2, P1)
# formula
return np.array( [vx1, vy1,
gs1x + g12x,
gs1y + g12y,
vx2, vy2,
gs2x + g21x,
gs2y + g21y] )
t = np.arange(0, 30, 0.01) # time
# solve ode vector eq
# 太陽座標
sx, sy = 0, 0
M = 10000
m1, m2 = 200, 100
r1, theta1 = 3.5, 220
r2, theta2 = 6, 135
x1, y1 = polarCoordinate(r1, theta1)
x2, y2 = polarCoordinate(r2, theta2)
# 圓形軌道速度
v1 = circleOrbitSpeed(M, r1)
v2 = circleOrbitSpeed(M, r2)
vx1, vy1 = tangentVector(v1, theta1)
vx2, vy2 = tangentVector(v2, theta2)
# 逃逸速度 C = 2^0.5 = 1.4142135623730951
C1 = 1.0
C2 = 1.3
track = odeint(ode, ( x1, y1, C1*vx1, C1*vy1,
x2, y2, C2*vx2, C2*vy2 ),
t,
args = ((sx, sy), M, m1, m2))
K1 = K_Energy(m1, track[:,2], track[:,3])
K2 = K_Energy(m2, track[:,6], track[:,7])
Rs1 = norm(track[:,0], track[:,1])
Rs2 = norm(track[:,4], track[:,5])
R12 = norm(track[:,0]-track[:,4], track[:,1]-track[:,5])
U12 = U_Energy(m1, m2, R12)
U1 = U_Energy(M, m1, Rs1) + U12
U2 = U_Energy(M, m2, Rs2) + U12
# 運行軌跡
plt.figure()
plt.plot(track[:,0], track[:,1], label='m1 = %s' %m1)
plt.plot(track[:,4], track[:,5], label='m2 = %s' %m2)
plt.plot(sx, sy, 'ko')
plt.title("sun mass M: %s" %M)
plt.xlabel('phase x')
plt.ylabel('phase y')
plt.legend()
# 動能與位能
plt.figure()
plt.plot(t, K1, label='K1')
plt.plot(t, K2, label='K2')
plt.plot(t, U1, label='U1')
plt.plot(t, U2, label='U2')
plt.legend()
# 能量守恆
plt.figure()
E1 = K1 + U1
E2 = K2 + U2
E = E1 + E2
plt.plot(t, E1, label='E1')
plt.plot(t, E2, label='E2')
plt.plot(t, E, label='E')
plt.legend()
plt.show()
--
※ 發信站: 批踢踢實業坊(ptt.cc)
◆ From: 111.251.155.111
推
03/27 21:21, , 1F
03/27 21:21, 1F
推
03/27 23:47, , 2F
03/27 23:47, 2F
※ 編輯: beatitude 來自: 111.251.152.54 (04/03 23:41)
討論串 (同標題文章)