Python で微分方程式の数値解法
Python
2017 10 25
1
x(t) = (x1(t), . . . , xn(t)) x(t0) = (x1(t0), . . . , xn(t0)) M
dx1 dt
=
v1(x, t),
...
x1(t0) = x10
dxn dt
=
vn(x, t),
xn(t0) = xn0
dx = v(x, t),
dt
x(t0) = x0
x t v Rn v : M Rn (vector filed)
t x
dx
x(t + t) - x(t)
= lim
dt t0
t
t
x(t + t) - x(t)
t
f (x, t)
t t x(t)
t = t0 x(t0) {xn}(n = 0, 1, 2, . . . ) x0 = x(t0) t
R3 Lorentz
dx dt = -x + y dy dt = rx - y + xz dz dt = -bz + xy
Lorenz (1963) = 10, r = 24.74, b = 8/3
1.1 Euler
Euler x0, x1, . . . , xn, . . .
xn+1 = xn + tf (xn, tn), tn = t0 + nt.
1
Python euler_orbit(x0, T, dt, f) Euler vectorfield v m [v1, . . . , vm]
t0 m 2 width x0 = [x01, . . . , x0m] t0 t0 + T dt N = T /dt m [x0, x1, . . . , xN ] 3 x0 x while orbit.append(list(x)) orbit
1 def euler_orbit(x0, T, dt, vectorfield):
2
width = len(x0)
3
x = x0
4
t=0
5
orbit = []
6
orbit.append(list(x))
7
while t > def func(x, y, *args):
>>
print a, b, args
>> print func('a', 'b', 'c', 'd', 'e', 'f')
('a', 'b', ('c', 'd', 'e', 'f'))
>> print func('a', 'b', 'c', 'd')
('a', 'b', ('c', 'd'))
>> print func('a', 'b', 'c')
2
('a', 'b', ('c',)) >> print func('a', 'b') ('a', 'b', ())
1.3 Euler Lorenz
Lorenz (25), (26), (27) Euler diff_euler.py
vector = [v1, v2, v3] Lorenz x0 = [0.1, 0.1, 0.1]
1 0.01 Euler 100 orbit
1 #!/usr/bin/env python
diff euler.py
2 # -*- coding: utf-8 -*-
3
4 # vector fields v1, v2, v3 of Lorenz eq as functions of t, x, y, z.
5 p, r, b = 10.0, 28.0, 8.0/3.0
6 def v1(t, x, y, z):
7
return(-p * x + p * y)
8 def v2(t, x, y, z):
9
return(-x * z + r * x - y)
10 def v3(t, x, y, z):
11
return(x * y - b * z)
12
13 def euler_orbit(x0, T, dt, vectorfield):
14
width = len(x0)
15
x = x0
16
t=0
17
orbit = []
18
while t >> import numpy as np 2 >>> a = [[1,2,3], [4,5,6], [7,8,9], [10,11,12]] 3 >>> a[1][2] 46 5 >>> a[:][2] 6 [7, 8, 9] 7 >>> na = np.array(a) 8 >>> na[:, 2] 9 array([ 3, 6, 9, 12]) 10 >>> na[1, :]
4
11 array([4, 5, 6])
7 na = np.array(a) ndarray 8 3 10 2 array ndarry A A[j, k] A[j][k]
NumPy orbit 3 nporbit = np.array(orbit) N ? 3 ndarray xorbit, yorbit, zorbit 16 matplotlib
diff euler.py 1 .... 2 orbit = euler_orbit(x0, T, dt, vector) 3 nporbit = np.array(orbit)
4
5 #xorbit = [] 6 #yorbit = [] 7 #zorbit = [] 8 #for x, y, z in orbit: 9 # xorbit.append(x) 10 # yorbit.append(y) 11 # zorbit.append(z)
12
13 # plot in 3-din space 14 fig = plt.figure() 15 ax = fig.gca(projection='3d') 16 ax.plot(nporbit[:, 0], nporbit[:, 1], nporbit[:, 2]) 17 .....
1.3 Euler NumPy ndarray diff_euler.py T= 50, dt = 0.01
1.4 Runge-Kutta
Euler Runge-Kutta dt
Python runge_kutta_orbit(x0, T, dt, vectorfield) Runge-Kutta vectorfield v m [v1, . . . , vm] t0 m 2 width x0 = [x01, . . . , x0m] t0 t0+T dt N = T /dt m [x0, x1, . . . , xN ] 3 x0 x while orbit.append(list(x)) orbit
Euler while tn = nt xn tn+1 = tn + dt xn+1 22,23 dt/2 k1, k2, k3, k4 4
1 def runge_kutta(x0, T, dt, vectorfield):
2
width = len(x0)
3
x = x0
4
t = 0.0
5
orbit = []
6
while t ................
................
In order to avoid copyright disputes, this page is only a partial summary.
To fulfill the demand for quickly locating and searching documents.
It is intelligent file search solution for home and business.
Related searches
- python replace multiple character in string
- python yield function
- yahoo finance python module
- philosophers song monty python lyrics
- monty python drinking song lyrics
- monty python the philosophers song
- python venomous bite
- ball python venomous
- python list of dictionaries append
- python thousands separator python 3
- python including other python files
- python get python file path