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

252 8 SolvingOrdinaryDifferentialEquations # Compute the time points where we want the solution N_t = int(round(T/dt)) time_points = np.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 plt.plot(t[-m:], u[-m:]) plt.hold(’on’) legends.append(solver.name()) plt.xlabel(’t’) # Plot exact solution too plt.plot(t[-m:], X_0*np.cos(omega*t)[-m:], ’k--’) legends.append(’exact’) plt.legend(legends, loc=’lower left’) plt.axis([t[-m], t[-1], -2*X_0, 2*X_0]) plt.title(’Simulation of {:d} periods with {:d} intervals per period’\ .format(number_of_periods, time_intervals_per_period)) plt.savefig(’tmp.pdf’); plt.savefig(’tmp.png’) plt.show() Anewfeature in thiscode is theability toplotonly the lastpperiods,whichallows us to perform long time simulations and watch the end results without a cluttered plot with too many periods. The syntax t[-m:] plots the last m elements in t (a negative index inPythonarrays/listscounts fromtheend). We may compareHeun’s method (i.e., the RK2 method) with the Euler-Cromer scheme: compare(odespy_methods=[odespy.Heun, odespy.EulerCromer], omega=2, X_0=2, number_of_periods=20, time_intervals_per_period=20) Figure 8.25 shows how Heun’s method (blue line) has considerable error in both amplitudeand phase alreadyafter 14–20periods (upper left), but using three times as many time steps makes the curves almost equal (upper right). However, after 194–200periods the errorshavegrown(lower left), but can be sufficiently reduced byhalving the time step (lower right). With all the methods in Odespy at hand, it is now easy to start exploring other methods, such as backward differences instead of the forward differences used in theForwardEuler scheme.Exercise8.22addresses thatproblem. Odespycontainsquitesophisticatedadaptivemethodswhere theuser is“guaran- teed”togetasolutionwithprescribedaccuracy.Thereisnomathematicalguarantee, but the error will for most cases not deviate significantly from the user’s tolerance that reflects the accuracy. A very popular method of this type is the Runge-Kutta- Fehlberg method, which runs a fourth-orderRunge-Kutta method and uses a fifth- order Runge-Kutta method to estimate the error so that Δt can be adjusted to keep the error below a tolerance. This method is also widely known as ode45, becausethat is thenameofthefunctionimplementingthemethodinMatlab.Wecan