In Runge Kutta Methods Why Midpont Method is accurate than Modified Euler Method?

522 Views Asked by At

Both Midpoints Method and Modified Euler Method is local truncation error

$ O(h^2) $

but Midpoint Method is accurate than Modified Euler Method

why?

1

There are 1 best solutions below

0
On

To get an empirical picture of the error sizes, take some non-linear test problem $$ y'+10\sin(y)=p'+10\sin(p)\iff y'=-10(\sin(y)-\sin(p))+p' $$ with $p(x)=\sin(x)$ as exact solution and apply the trapezoidal and midpoint methods in the explicit and implicit variants over a couple of step sizes. To get comparable quantities, compute the difference to the exact solution and divide by the expected size $h^2$. This gives the plots

enter image description here

As one can see, the errors for the midpoint methods have about the same magnitude. The explicit trapezoidal method has an error of the double size, while the implicit trapezoidal method has one tenth of the error compared to the midpoint methods.

The implementation for the methods is

def exp_mid(f,x,y,h): k=h*f(x,y); return y+h*f(x+0.5*h,y+0.5*k);

def exp_trap(f,x,y,h): k=h*f(x,y); return y+0.5*(k+h*f(x+h,y+k));

def imp_mid(f,x,y,h):
    k=h*f(x,y);
    for _ in range(6): k=h*f(x+0.5*h,y+0.5*k);
    return y+k;

def imp_trap(f,x,y,h):
    k1=h*f(x,y); k2=k1;
    for _ in range(6): k2=h*f(x+h,y+0.5*(k1+k2));
    return y+0.5*(k1+k2);

and the test data

def p(x): return np.sin(x)
def dp(x): return np.cos(x)

def ode(x,y): return -10*(np.sin(y)-np.sin(p(x)))+dp(x)

x0,xf=1,8;
steps=[0.02, 0.05, 0.1]

The plot is then produced via

fig, ax = plt.subplots(4,1,figsize=(8,4*2))
methods=[exp_mid, imp_mid, exp_trap, imp_trap]
names=["expl. midpoint", "impl. midpoint", "expl. trapez", "impl. trapez"]
for m,method in enumerate(methods):
    for h in steps:
        x = np.arange(x0,xf+h,h);
        y = [p(x0)]
        for k in range(len(x)-1):
            y.append(method(ode,x[k],y[k],h))
        y = np.asarray(y)
        ax[m].plot(x,(y-p(x))/h**2,'-o',ms=2+10*h, label="h=%.4f"%h)
    ax[m].legend(); ax[m].grid(); ax[m].set_title(names[m])
plt.tight_layout(); plt.show()