We'll solve the same problem again, using the numerical ODE solver from the GSL.
def f_1(t, y, params):
return [y[1], -y[0] - params[0] * y[1] * (y[0]**2 - 1.0)]
def j_1(t, y, params):
return [ [0.0, 1.0], [-2.0 * params[0] * y[0] * y[1] - 1.0,
-params[0] * (y[0] * y[0] - 1.0)] ]
T = ode_solver()
T.algorithm = "rk8pd"
T.function = f_1
T.jacobian = j_1
T.ode_solve(y_0=[2, 0], t_span=[0, 20], params=[1.0],
num_points=1000)
interpolator = T.interpolate_solution()
plot(interpolator, (0, 20), axes_labels=('t','y'), figsize=(4,3))The output is shown in the following screenshot:
This solver has a very different interface from the one in the previous example. The ordinary differential equations are defined by a Python function called f_1. The first argument to this function is the independent variable, the second is a list of dependent variables, and the third is a list of parameters. y[0] represents...