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

250 8 SolvingOrdinaryDifferentialEquations Odespy:ExamplewithOscillations UsingOdespytosolveoscillationODEslike u′′+ω2u= 0, reformulatedas a systemu′ = v andv′ =−ω2u, can be done with the followingcode: import odespy import numpy as np import matplotlib.pyplot as plt # Define the ODE system # u’ = v # v’ = -omega**2*u def f(sol, t, omega=2): u, v = sol return [v, -omega**2*u] # Set and compute problem dependent parameters omega = 2 X_0 = 1 number_of_periods = 40 time_steps_per_period = 20 P = 2*np.pi/omega # length of one period dt = P/time_steps_per_period # time step T = number_of_periods*P # final simulation time # Create Odespy solver object odespy_method = odespy.RK2 solver = odespy_method(f, f_args=[omega]) # The initial condition for the system is collected in a list solver.set_initial_condition([X_0, 0]) # Compute the desired time points where we want the solution N_t = int(round(T/dt)) # no of time intervals time_points = np.linspace(0, T, N_t+1) # Solve the ODE problem sol, t = solver.solve(time_points) # Note: sol contains both displacement and velocity # Extract original variables u = sol[:,0] v = sol[:,1] plt.plot(t, u, t, v) # ...for a quick check on u and v plt.show() After specifying the number of periods to simulate, as well as the number of time stepsperperiod,wecompute the time step (dt) andsimulationend time(T). Thetwostatementsu = sol[:,0]andv = sol[:,1]are important,sinceour twofunctionsuandv in theODEsystemarepackedtogether inonearrayinside the Odespysolver (thesolutionof the ODE system is returnedfromsolver.solveas a two-dimensionalarraywhere thefirstcolumn(sol[:,0])storesuandthesecond (sol[:,1])storesv).