Problem About the Order of Convergence of Numerical Method

Question

Sometimes when I did convergence test for my numerical solver, I can't get the exact order of convergence even I am sure that my code is correct. I would like to ask if anyone else has experienced this and knows what the cause is?

I put a very simple example here about the Euler Method. I want to approximate $x$ with $$ x'=\sqrt{x}+x $$ And I use the Euler Method such that $$ x(t_{0}+\Delta t) = x(t_{0})+x'(t_{0})\Delta t $$

I also attach my python code here:

delta_t_list = [0.001, 0.005, 0.01, 0.02]
error_list = []
for delta_t in delta_t_list:
    x_list = [1]
    max_iter_times = int(3/delta_t)
    t = np.linspace(0, 3, max_iter_times+1)
    x= (2*np.e**(t/2)-1)**2
    for i in range(max_iter_times):
        x0 = x_list[-1]
        x1 = x0+(np.sqrt(x0)+x0)*delta_t
        x_list.append(x1)
    error = np.linalg.norm(x-x_list, ord=2)
    error_list.append(error)

plt.plot(np.log(delta_t_list), np.log(error_list))
linregress(np.log(delta_t_list), np.log(error_list))

The result I got is

Answer 1

The norm you compute is just the Euclidean norm of the vector. As the sample number increases with decreasing step size, the norm increases with the square root of the length. You have to change to

    error = np.linalg.norm(x-x_list, ord=2)*delta_t**0.5

You might also need to account and correct for the possibility that the step size of the t array is not delta_t due to floating point errors in computing 3/delta_t. The safe way would be to set either

    t = np.arange(0,3,delta_t) # or 3+0.5*delta_t
    max_iter_times = len(t)-1

or to hit the end point correctly

    max_iter_times = int(round(3/delta_t))
    t, delta_t = np.linspace(0, 3, max_iter_times+1, retstep=True)