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

134 4 SolvingOrdinaryDifferentialEquations a = 2 b = 1 solver = method(f, f_args=[a, b]) This is a good feature because problemparametersmust otherwise be global vari- ables–nowtheycanbearguments inour right-handside function inanaturalway. Exercise4.16asksyoutomakeacomplete implementationof thisproblemandplot thesolution. UsingOdespy to solve oscillationODEs likeu00 C!2u D 0, reformulated as a systemu0 D v andv0 D !2u, is done as follows. We specify a given number of time stepsperperiodandcompute the associated time steps andend timeof the simulation (T),givenanumberofperiods to simulate: import odespy # 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_intervals_per_period = 20 from numpy import pi, linspace, cos P = 2*pi/omega # length of one period dt = P/time_intervals_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 = 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] The last two statements are important sinceour two functionsu andv in theODE system are packed together in one array inside the Odespy solver. The solution