Page - 135 - in Programming for Computations – Python - A Gentle Introduction to Numerical Simulations with Python

4.3 OscillatingOne-DimensionalSystems 135 of theODE system is returned as a two-dimensional arraywhere the first column (sol[:,0]) stores u and the second (sol[:,1]) stores v. Plotting u and v is amatterof runningplot(t, u, t, v). Remark In the right-hand side functionwewrite f(sol, t, omega) instead of f(u, t, omega) to indicate that the solution sent tof is a solution at timetwhere thevaluesofuandv arepacked together:sol = [u, v].Wemightwelluseu asargument: 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 solutionat timet for thefirst componentof theODEsystem. To switch to another numericalmethod, just substituteRK2by the propername of the desiredmethod. Typing pydoc odespy in the terminal windowbrings up a list of all the implementedmethods. Thisvery simplewayofchoosingamethod suggests anobvious extension of the code above: we candefine a list ofmethods, run all methods, and compare theiru curves in a plot. As Odespy also contains theEuler-Cromerscheme,we rewrite the systemwithv0 D !2uas thefirstODE andu0 Dv as the secondODE,because this is the standardchoicewhenusing 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. Thecodeosc_odespy.pycontains thedetails: def compare(odespy_methods, omega, X_0, number_of_periods, time_intervals_per_period=20): from numpy import pi, linspace, cos P = 2*pi/omega # length of one period dt = P/time_intervals_per_period T = number_of_periods*P