Error upper bound using Euler’s Method

euler's methodinitial-value-problemsnumerical methodsordinary differential equations

Determine an upper bound on the error made using Euler's method with step size $h$ to find an approximate value of the solution to the initial-value problem:

$\frac{dy}{dt} = t – y^4$, $y(0) = 0$

at any point $t$ in the interval $[0, 1]$.

Knowing that $f(t, y) = \frac{dy}{dt} = t – y^4$, I calculated $\frac{\partial f}{\partial y} = -4y^3$. Unsure where to go from here.

Best Answer

Let's look at the half axis $y=0$, $t>0$. There the right side is $f(t,0)=t>0$ so that no solution may cross from the upper to the lower quadrant.

Close to zero one gets $y(t)=\frac12t^2+O(t^9)$ so that the solution will indeed enter the upper quadrant from the start. Furthermore, from $y'(t)\le t$ we get $y(t)\le\frac12t^2$, so that we also know an upper bound for the solution. We can restrict the region for the estimates of the Euler method to $(t,x)\in[0,1]\times[0,1]$, or, if you want to be cautious, $(t,x)\in[0,1]\times[-1,1]$.

On that region, $$|f(t,y)|\le 1=M_1$$ is a bound for the first derivative of any solution, and $$|f_t+f_yf|=|1-4y^3(t-y^4)|\le 5=M_2$$ a bound for the second derivative. Also, for the $y$-Lipschitz constant one gets similarly $$|f_y|=|-4y^3|\le 4=L. $$

Now insert into the error estimate $$ e(t,h)\le \frac{M_2}{2L}(e^{Lt}-1)h=\frac{5}{8}(e^{4t}-1)h. $$

Now let's see how that bound stands up to the actual error of the numerical method

def EulerIntegrate(f,y0,t):
    y = np.asarray(len(t)*[y0], dtype=float);
    for k in range(len(t)-1):
        y[k+1] = y[k]+(t[k+1]-t[k])*f(t[k],y[k])
    return y

def f(t,y): return t-y**4

t = np.linspace(0,1,64*15+1)

N = [30*2**k for k in range(4)]
yEuler = [ EulerIntegrate(f,0,t[::n]) for n in N];

yGold = odeint(f,0,t,tfirst=True, atol=1e-12, rtol=1e-13)[:,0]

fig,ax = plt.subplots(1,2,figsize=(12,5))
for k in range(len(N)):
    h = t[N[k]]-t[0];
    ax[0].plot(t[::N[k]], yEuler[k],'-o',label="h=%.4f"%(1./N[k]))
    ax[1].plot(t[::N[k]], (yEuler[k]-yGold[::N[k]])/h,'-o',label="h=%.4f"%(1./N[k]))

gives the plots

enter image description here

where the second plot shows the error profile, the estimated leading coefficient $c(t)$ of the global error $e(t,h)=c(t)h+O(h^2)$ over time. The rapidly falling gray line is the error bound, safely below the actual error. The left plot of the actual solutions against the backdrop of a much more precise numerical solution clearly shows the linear convergence of the Euler method.

Related Question