课程概要：
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时差距就会很明显，但第二个图与视频中绘出的图有很大区别
