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

136 4 SolvingOrdinaryDifferentialEquations # 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]) # Compute the time points where we want the solution dt = float(dt) # avoid integer division N_t = int(round(T/dt)) time_points = linspace(0, N_t*dt, N_t+1) legends = [] for solver in solvers: sol, t = solver.solve(time_points) v = sol[:,0] u = sol[:,1] # Plot only the last p periods p = 6 m = p*time_intervals_per_period # no time steps to plot plot(t[-m:], u[-m:]) hold(’on’) legends.append(solver.name()) xlabel(’t’) # Plot exact solution too plot(t[-m:], X_0*cos(omega*t)[-m:], ’k--’) legends.append(’exact’) legend(legends, loc=’lower left’) axis([t[-m], t[-1], -2*X_0, 2*X_0]) title(’Simulation of %d periods with %d intervals per period’ % (number_of_periods, time_intervals_per_period)) savefig(’tmp.pdf’); savefig(’tmp.png’) show() Anewfeature in thiscode is theability toplotonly the lastpperiods,whichallows us to perform long time simulations andwatch the end resultswithout a cluttered plot with too many periods. The syntax t[-m:] plots the last m elements in t (anegative index inPythonarrays/lists counts fromtheend). Wemay compare Heun’s method (or equivalently the RK2method) with the Euler-Cromerscheme: compare(odespy_methods=[odespy.Heun, odespy.EulerCromer], omega=2, X_0=2, number_of_periods=20, time_intervals_per_period=20) Figure4.22showshowHeun’smethod (theblue linewith small disks) has consid- erable error in both amplitude and phase already after 14–20 periods (upper left), but using three times as many time steps makes the curves almost equal (upper