Page - 251 - in Programming for Computations – Python - A Gentle Introduction to Numerical Simulations with Python 3.6, Volume Second Edition

8.4 Oscillating1DSystems:ASecondOrderODE 251 Remark In the right-handsidefunctionwewritef(sol, t, omega) insteadoff(u, t, omega) to indicate that thesolutionsent tof is a solutionat timetwhere thevaluesofuandvarepackedtogether:sol = [u, v].Wemightwelluse uasargument: def f(u, t, omega=2): u, v = u return [v, -omega**2*u] This just means that we redefine the nameu inside the function to mean the solutionat timet for thefirst componentof theODEsystem. To switch to another numerical method, just substituteRK2by the proper name of the desired method. Typingpydoc 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 theiru curves in a plot. As Odespy also contains the Euler-Cromer scheme, we rewrite the system with v′ = −ω2u as the first ODE andu′ = v as the second ODE, because this is the standard choice when using the Euler-Cromermethod(also inOdespy): 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 functioncompare inosc_odespy.pycontains 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])