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

272 8 SolvingOrdinaryDifferentialEquations With (8.87),we havea simplewayofcomputing theerrorEi, letting Ei =‖en‖l2 . (8.88) We are nowin position to computeconvergencerates, andwrite correspondingtest functions,also forODEsolvers. 8.5.3 TestFunction:ConvergenceRatesfortheFESolver To illustrate, we write a simple test function for ode_FE that we implemented previously(Sect. 8.2.5).Applying the solver to a populationgrowthmodel, the test functioncouldbewritten: def test_convergence_rates_ode_FE(number_of_experiments): """ Test that the convergence rate with the ode_FE solver is 1. Use population growth model as test case. """ U_0=100 # initial value T=20 # total time span dt = 2.0 # initial time step expected_rate_FE = 1.0 def f(u, t): """Population growth, u’ = a*u, with a = 0.1 here.""" return 0.1*u def u_exact(t): return 100*np.exp(0.1*t) dt_values = [] E_values = [] for i in range(number_of_experiments): u, t = ode_FE(f=f, U_0=U_0, dt=dt, T=T) u_e = u_exact(t) # get exact solution at time mesh points (in t) E = np.sqrt(dt*np.sum((u_e-u)**2)) # ...discrete L^2 norm dt_values.append(dt) E_values.append(E) dt = dt/2 # Halving time step for next solve r = [np.log(E_values[i]/E_values[i-1])/ np.log(dt_values[i]/dt_values[i-1]) for i in range(1, number_of_experiments, 1)] #print(r) # Accept rate to 1 decimal place tol = 0.1 assert abs(r[-1] - expected_rate_FE) < tol return When test_convergence_rates_ode_FE is called, the for loop will execute ode_FE the number of times specified by the input parameter number_of_experiments. Each execution of ode_FE happens with half the time step of the previous execution. Errors (E) and time steps (dt) are stored in