import numpy as np import matplotlib.pyplot as plt
然后,我们从将用于这些方程的变量开始。
sigma =10.0 rho =28.0 beta =8/3 t =0 tf =40 h =0.01
设置好变量后,我们现在继续讨论导数函数。
defderivative(r,t): x = r[0] y = r[1] z = r[2]return np.array([sigma * (y - x), x * (rho - z) - y, (x * y) - (beta * z)])
快速回顾一下导数函数。 r 变量用于在数组内设置初始条件。 我们为每个变量分配数组的值,在本例中为三个,然后函数将返回开始时所示的微分方程,但现在具有初始条件。 每次循环经过一次迭代,函数将继续更新新值,这些值将存储在以下数组中。
time = np.array([]) x = np.array([]) y = np.array([]) z = np.array([]) r = np.array([1.0, 1.0, 1.0])
现在,我们的主要部分设置完毕,我们只需继续将所有内容应用到 RK4 循环内,我们最终会得到
while (t <= tf ): time = np.append(time, t) z = np.append(z, r[2]) y = np.append(y, r[1]) x = np.append(x, r[0]) k1 = h*derivative(r,t) k2 = h*derivative(r+k1/2,t+h/2) k3 = h*derivative(r+k2/2,t+h/2) k4 = h*derivative(r+k3,t+h) r += (k1+2*k2+2*k3+k4)/6 t = t + h