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.

Google Online Preview   Download