python科学计算最佳实践scipy指南_10-Python 科学计算_scipy 篇

课程概要：

1、非线性方程组

2、数值积分

3、常微分方程组

1、非线性方程组

'''

5 * x1 - 25 = 0

5 * x0 * x0 - x1 * x2 = 0

x2 * x0 -27 = 0

scipy.optimize fsolve(func, x)

# 所使用的scipy中的库optimize以及方法fsolve

# func 是自己构造的函数

'''

from scipy.optimize import fsolve

def func(x):

x0, x1, x2 = x.tolist()

return [

5 * x1 - 25,

5 * x0 * x0 - x1 * x2,

x2 * x0 -27

]

r = fsolve(func,[1,1,1])

print r

>>>

[ 3. 5. 9.]

'''

5 * x1 + 3 = 0

5 * x0 * x0 - 2sin(x1 * x2) = 0

x1 * x2 -1.5 = 0

'''

from scipy.optimize import fsolve

from math import sin

def func(x):

x0, x1, x2 = x.tolist()

return [

5 * x1 + 3,

5 * x0 * x0 - 2 * sin(x1 * x2) ,

x1 * x2 -1.5

]

r = fsolve(func,[1,1,1])

print r

>>>

[-0.63166288 -0.6 -2.5 ]

'''

5 * x1 + 3 = 0

5 * x0 * x0 - 2sin(x1 * x2) = 0

x1 * x2 -1.5 = 0

'''

from scipy.optimize import fsolve

from math import sin, cos

def func(x):

x0, x1, x2 = x.tolist()

return [

5 * x1 + 3,

5 * x0 * x0 - 2 * sin(x1 * x2) ,

x1 * x2 -1.5

]

def j(x):

x0, x1, x2 = x.tolist()

return [

[0, 5, 0], # 分别对x0, x1, x2 求偏导

[10* x0, -2 * x2 * cos(x1 * x2), -2 * x1 * cos(x1 * x2)],

[0, x2, x1]

]

r = fsolve(func, [1,1,1], fprime=j) # 使用求导矩阵传入时可以提高4倍效率

print r

>>>

[-0.63166288 -0.6 -2.5]

2、数值积分(定积分)

# 求半圆面积(定义求)

'''

pi * 1 *1 / 2

x * x + y * y = 1

y = (1 - x * x) ** 0.5

'''

import numpy as np

def func(x):

return (1 - x * x) ** 0.5

N = 10000

x = np.linspace(-1, 1, N)

dx = 2.0 / N

y = func(x)

m = dx * np.sum(y[:-1] + y[1:])

print m / 2 # m是整圆的面积

>>>

1.570637584

# 使用scipy的库进行计算

import numpy as np

from scipy import integrate

def func(x):

return (1 - x * x) ** 0.5 # x^2 + y^2 = 1

p, err = integrate.quad(func, -1, 1) # x 给定的区间

print p

>>>

1.57079632679

# 使用scipy的库进行计算

import numpy as np

from scipy import integrate

def func(x):

return (4 - x * x) ** 0.5 # x^2 + y^2 = 4

p, err = integrate.quad(func, -2, 2)

print p

>>>

6.28318530718

# dblquad：双重积分 trlquad：三重积分

'''

x * x+ y * y + z * z =1

z = (1 - x * x - y * y)**0.5

'''

from scipy import integrate

def func(x, y):

return (1 - x * x - y * y)**0.5

def func2(x):

return (1 - x * x)**0.5

m = integrate.dblquad(func, -1, 1, # x 的区间

lambda x: - func2(x),

lambda x: func2(x)) # y 的积分区间

print m

>>>

(2.0943951023931984, 2.3252456653390912e-14)

3、常微分方程组

dx/dt=σ(y-x)

dy/dt=x(ρ-z)-y

dz/dt=xy-βz

# σ = p , ρ = r , β= b

from scipy.integrate import odeint

import numpy as np

def lorenz(w, t, p, r, b):

x, y, z = w

return np.array([p * (y-x), x * (r-z) - y, x * y - b * z])

t = np.arange(0, 30, 0.01)

track1 = odeint(lorenz, (0.0, 1.00, 0.0), t, args=(10.0, 28.0, 3.0)) # 初始值，即w

## track2 = odeint(lorenz, (0.0, 1.01, 0.0), t, args=(10.0, 28.0, 3.0)) # args 为(p, r, b)

from mpl_toolkits.mplot3d import Axes3D

import matplotlib.pyplot as plt

fig = plt.figure()

ax = Axes3D(fig)

ax.plot(track1[:, 0], track1[:, 1], track1[:, 2])

plt.show()

# 改变轨迹值

ax.plot(track2[:, 0], track2[:, 1], track2[:, 2])

# 只有细微的差别，当改为1.10时差距就会很明显，但第二个图与视频中绘出的图有很大区别