def convergence_rates(m, num_periods=8):
"""
Return m-1 empirical estimates of the convergence rate
based on m simulations, where the time step is halved
for each simulation.
"""
w = 0.35; I = 0.3
dt = 2*pi/w/30 # 30 time step per period 2*pi/w
T = 2*pi/w*num_periods
dt_values = []
E_values = []
for i in range(m):
u, t = solver(I, w, dt, T)
u_e = exact_solution(t, I, w)
E = sqrt(dt*sum((u_e-u)**2))
dt_values.append(dt)
E_values.append(E)
dt = dt/2
r = [log(E_values[i-1]/E_values[i])/
log(dt_values[i-1]/dt_values[i])
for i in range(1, m, 1)]
return r
Result: r
contains values equal to 2.00 - as expected!