8.4 Oscillating 1D Systems: A Second Order ODE
251
Remark
In the right-hand side function we write f(sol, t, omega) instead of f(u,
t, omega) to indicate that the solution sent to f is a solution at time t where
the values of u and v are packed together: sol = [u, v]. We might well use
u as argument:
def f(u, t, omega=2):
u, v = u
return [v, -omega**2*u]
This just means that we redefine the name u inside the function to mean the
solution at time t for the first component of the ODE system.
To switch to another numerical method, just substitute RK2 by the proper name
of the desired method. Typing pydoc odespy in the terminal window brings up a
list of all the implemented methods. This very simple way of choosing a method
suggests an obvious extension of the code above: we can define a list of methods,
run all methods, and compare their u curves in a plot. As Odespy also contains the
Euler-Cromer scheme, we rewrite the system with v = −ω 2 u as the first ODE
and u = v as the second ODE, because this is the standard choice when using the
Euler-Cromer method (also in Odespy):
def f(u, t, omega=2):
v, u = u
return [-omega**2*u, v]
This change of equations also affects the initial condition: the first component
is zero and second is X_0, so we need to pass the list [0, X_0] to
solver.set_initial_condition.
The function compare in osc_odespy.py contains the details:
def compare(odespy_methods,
omega,
X_0,
number_of_periods,
time_intervals_per_period=20):
P = 2*np.pi/omega
# length of one period
dt = P/time_intervals_per_period
T = number_of_periods*P
# If odespy_methods is not a list, but just the name of
# a single Odespy solver, we wrap that name in a list
# so we always have odespy_methods as a list
if type(odespy_methods) != type([]):
odespy_methods = [odespy_methods]
# Make a list of solver objects
solvers = [method(f, f_args=[omega]) for method in
odespy_methods]
for solver in solvers:
solver.set_initial_condition([0, X_0])
251
Remark
In the right-hand side function we write f(sol, t, omega) instead of f(u,
t, omega) to indicate that the solution sent to f is a solution at time t where
the values of u and v are packed together: sol = [u, v]. We might well use
u as argument:
def f(u, t, omega=2):
u, v = u
return [v, -omega**2*u]
This just means that we redefine the name u inside the function to mean the
solution at time t for the first component of the ODE system.
To switch to another numerical method, just substitute RK2 by the proper name
of the desired method. Typing pydoc odespy in the terminal window brings up a
list of all the implemented methods. This very simple way of choosing a method
suggests an obvious extension of the code above: we can define a list of methods,
run all methods, and compare their u curves in a plot. As Odespy also contains the
Euler-Cromer scheme, we rewrite the system with v = −ω 2 u as the first ODE
and u = v as the second ODE, because this is the standard choice when using the
Euler-Cromer method (also in Odespy):
def f(u, t, omega=2):
v, u = u
return [-omega**2*u, v]
This change of equations also affects the initial condition: the first component
is zero and second is X_0, so we need to pass the list [0, X_0] to
solver.set_initial_condition.
The function compare in osc_odespy.py contains the details:
def compare(odespy_methods,
omega,
X_0,
number_of_periods,
time_intervals_per_period=20):
P = 2*np.pi/omega
# length of one period
dt = P/time_intervals_per_period
T = number_of_periods*P
# If odespy_methods is not a list, but just the name of
# a single Odespy solver, we wrap that name in a list
# so we always have odespy_methods as a list
if type(odespy_methods) != type([]):
odespy_methods = [odespy_methods]
# Make a list of solver objects
solvers = [method(f, f_args=[omega]) for method in
odespy_methods]
for solver in solvers:
solver.set_initial_condition([0, X_0])