常微分方程解法（2）例题

用SciPyのscipy.integrate.odeint解联立微分方程式，以时间t（tの范围在0~几秒就可以，例：t=0~2秒）和解x(t), y(t), z(t)作图。

作为t=0的初期条件，x(t=0) = -10, y (t=0) = 0, z (t=0) = 35.0。还有，关于系数a,b,c，可以尝试a = 40, b = 5, c = 35 的情况和a = 40, b = 10, c = 35的情况。还有时间分割Δt要适当取小一点的值。另外，Δt如果取的过小，计算量会就变得很大。根据实际实行后的错误进行调整吧。

下面是解决的python程序：

import numpy as np
import scipy.integrate as sciin
import matplotlib.pyplot as plt

#把与自变量t相关的导函数的函数名字放在F里
def f(F, t,params):    
    x,y,z = F       
    f_values = [a*(y-x),(c-a)*x-x*z+c*y,x*y-b*z]   #分别写x,y,z的导数等于的右边的式子 
    return f_values


# 放系数
a = 40           
b= 5                
c = 35       
#把上面三个系数放到parameters里
parameters = [a,b,c]

#设定x,y,z的初始值
x0 = -10      
y0 = 0.0  
z0 = 35  
#　把初始值放到Y0里
Y0 = [x0,y0,z0]

# 开始点，结束点，间隔的设置
tStart = 0.0        
tStop  = 2   
tInc   = 0.01       #间隔
#把上面归纳到t里
t = np.arange(tStart, tStop, tInc)

# sciin.odeint解ODE
solution = sciin.odeint(f, Y0, t, args=(parameters,))

# 作图
plt.figure(figsize=(9.5, 6.5))
plt.plot(t, solution[:, 0], color='black')
plt.plot(t, solution[:, 1], color='green')
plt.plot(t, solution[:, 2], color='red')
plt.xlabel('time, t' , fontsize=14)
plt.ylabel('theta(t)', fontsize=14)
plt.show()

 我简直是天才。