5.scipy Integrate

1import pylab as pl
2import numpy as np
3from scipy import integrate
4from scipy.integrate import odeint
5import matplotlib as mpl
6mpl.rcParams['font.sans-serif'] = ['SimHei']

1def half_circle(x):
2    return (1-x**2)**0.5

1N = 10000
2x = np.linspace(-1, 1, N)
3dx = x[1] - x[0]
4y = half_circle(x)
52 * dx * np.sum(y) # 面积的两倍 

1np.trapz(y, x) * 2 # 面积的两倍

1from scipy import integrate
2pi_half, err = integrate.quad(half_circle, -1, 1)
3pi_half * 2

1def half_sphere(x, y):
2    return (1-x**2-y**2)**0.5

1volume, error = integrate.dblquad(half_sphere, -1, 1, 
2        lambda x:-half_circle(x), 
3        lambda x:half_circle(x))
4
5print (volume, error, np.pi*4/3/2)

 1#%fig=洛伦茨吸引子：微小的初值差别也会显著地影响运动轨迹
 2from scipy.integrate import odeint 
 3import numpy as np 
 4
 5def lorenz(w, t, p, r, b): #❶
 6    # 给出位置矢量w，和三个参数p, r, b计算出
 7    # dx/dt, dy/dt, dz/dt的值
 8    x, y, z = w.tolist()
 9    # 直接与lorenz的计算公式对应 
10    return p*(y-x), x*(r-z)-y, x*y-b*z
11
12t = np.arange(0, 30, 0.02) # 创建时间点 
13# 调用ode对lorenz进行求解, 用两个不同的初始值 
14track1 = odeint(lorenz, (0.0, 1.00, 0.0), t, args=(10.0, 28.0, 3.0)) #❷
15track2 = odeint(lorenz, (0.0, 1.01, 0.0), t, args=(10.0, 28.0, 3.0)) #❸
16#%hide
17from mpl_toolkits.mplot3d import Axes3D
18fig = pl.figure()
19ax = Axes3D(fig)
20ax.plot(track1[:,0], track1[:,1], track1[:,2], lw=1)
21ax.plot(track2[:,0], track2[:,1], track2[:,2], lw=1);

1def mass_spring_damper(xu, t, m, k, b, F):
2    x, u = xu.tolist()
3    dx = u
4    du = (F - k*x - b*u)/m
5    return dx, du 

 1#%fig=滑块的速度和位移曲线
 2m, b, k, F = 1.0, 10.0, 20.0, 1.0
 3init_status = 0.0, 0.0
 4args = m, k, b, F
 5t = np.arange(0, 2, 0.01)
 6result = odeint(mass_spring_damper, init_status, t, args)
 7#%hide
 8fig, (ax1, ax2) = pl.subplots(2, 1)
 9ax1.plot(t, result[:, 0], label=u"位移")
10ax1.legend()
11ax2.plot(t, result[:, 1], label=u"速度")
12ax2.legend();

 1from scipy.integrate import ode
 2
 3class MassSpringDamper(object): #❶
 4    
 5    def __init__(self, m, k, b, F):
 6        self.m, self.k, self.b, self.F = m, k, b, F
 7        
 8    def f(self, t, xu):
 9        x, u = xu.tolist()
10        dx = u
11        du = (self.F - self.k*x - self.b*u)/self.m
12        return [dx, du] 
13
14system = MassSpringDamper(m=m, k=k, b=b, F=F)
15init_status = 0.0, 0.0
16dt = 0.01
17
18r = ode(system.f) #❷
19r.set_integrator('vode', method='bdf')
20r.set_initial_value(init_status, 0)
21
22t = []
23result2 = [init_status]
24while r.successful() and r.t + dt < 2: #❸
25    r.integrate(r.t + dt)
26    t.append(r.t)
27    result2.append(r.y)
28    
29result2 = np.array(result2)
30np.allclose(result, result2)

 1class PID(object):
 2    
 3    def __init__(self, kp, ki, kd, dt):
 4        self.kp, self.ki, self.kd, self.dt = kp, ki, kd, dt
 5        self.last_error = None
 6        self.status = 0.0
 7        
 8    def update(self, error):
 9        p = self.kp * error
10        i = self.ki * self.status
11        if self.last_error is None:
12            d = 0.0
13        else:
14            d = self.kd * (error - self.last_error) / self.dt
15        self.status += error * self.dt
16        self.last_error = error
17        return p + i + d

 1#%fig=使用PID控制器让滑块停在位移为1.0处
 2def pid_control_system(kp, ki, kd, dt, target=1.0):
 3    system = MassSpringDamper(m=m, k=k, b=b, F=0.0)
 4    pid = PID(kp, ki, kd, dt)
 5    init_status = 0.0, 0.0
 6
 7    r = ode(system.f)
 8    r.set_integrator('vode', method='bdf')
 9    r.set_initial_value(init_status, 0)
10
11    t = [0]
12    result = [init_status]
13    F_arr = [0]
14
15    while r.successful() and r.t + dt < 2.0:
16        r.integrate(r.t + dt)
17        err = target - r.y[0]  #❶
18        F = pid.update(err)  #❷
19        system.F = F  #❸
20        t.append(r.t)
21        result.append(r.y)
22        F_arr.append(F)
23
24    result = np.array(result)
25    t = np.array(t)
26    F_arr = np.array(F_arr)
27    return t, F_arr, result
28
29
30t, F_arr, result = pid_control_system(50.0, 100.0, 10.0, 0.001)
31print(u"控制力的终值:", F_arr[-1])
32#%hide
33fig, (ax1, ax2, ax3) = pl.subplots(3, 1, figsize=(6, 6))
34ax1.plot(t, result[:, 0], label=u"位移")
35ax1.legend(loc="best")
36ax2.plot(t, result[:, 1], label=u"速度")
37ax2.legend(loc="best")
38ax3.plot(t, F_arr, label=u"控制力")
39ax3.legend(loc="best")

 1%%time
 2from scipy import optimize
 3
 4
 5def eval_func(k):
 6    kp, ki, kd = k
 7    t, F_arr, result = pid_control_system(kp, ki, kd, 0.01)
 8    return np.sum(np.abs(result[:, 0] - 1.0))
 9
10
11kwargs = {
12    "method": "L-BFGS-B",
13    "bounds": [(10, 200), (10, 100), (1, 100)],
14    "options": {
15        "approx_grad": True
16    }
17}
18
19opt_k = optimize.basinhopping(
20    eval_func, (10, 10, 10), niter=10, minimizer_kwargs=kwargs)
21print(opt_k.x)

 1#%fig=优化PID的参数降低控制响应时间
 2kp, ki, kd = opt_k.x
 3t, F_arr, result = pid_control_system(kp, ki, kd, 0.01)
 4idx = np.argmin(np.abs(t - 0.5))
 5x, u = result[idx]
 6print ("t={}, x={:g}, u={:g}".format(t[idx], x, u))
 7#%hide
 8fig, (ax1, ax2, ax3) = pl.subplots(3, 1, figsize=(6, 6))
 9ax1.plot(t, result[:, 0], label=u"位移")
10ax1.legend(loc="best")
11ax2.plot(t, result[:, 1], label=u"速度")
12ax2.legend(loc="best")
13ax3.plot(t, F_arr, label=u"控制力")
14ax3.legend(loc="best");