[Math] Euler and Taylor methods and difference between Heun’s and Trapezoidal methods

calculusordinary differential equations

Heun (which is the improved Euler's method), 2nd order Taylor and Trapezoidal methods

I am a bit confused because I read two different things about the Heun's method and Trapezoidal method.

I cannot understand the difference between them, as I saw that in some places the Hen's method is the same as the explicit Trapezoidal method and in other place I saw that they are just very similar.

I am trying to find the explicit formula of the Trapezoidal method for $Y_{n+1}$
in terms of $t_n,t_{n+1},Y_n$ that uses the RHS of a certain function.

I have the function: $f(t,y)= \frac{8 \cos(4t)+y}{4}$

I found that the formula for the explicit Trapezoidal method is:

$Y_{n+1}=Y_n+\frac{h}{2}f(t_n,Y_n)+ \frac{h}{2}f(t_{n+1},Y_n+hf(t_n,Y_n)$

Because we use Euler's method to replace $Y_{n+1}$ on the RHS by:$ Y_{n+1}=Y_n+hf(t_n,Y_n) $ ( so, by doing this we will obtain the above formula)

Applying this formula in my function the answer will be :

$Y_{n+1}= Y_n+\frac{h}{2}[(\frac{8cos(4t_n)+Y_n}{4})+\frac{h}{2}(\frac{8cos(4t_{n+1})+Y_n}{4}+h(\frac{8cos(4t_n+Y_n}{4})]$

Now I need to know the formula for the Heun's method giving an explicit formula for $Y_{n+1}$ again in terms of $t_n,t_{n+1},Y_n$ that uses the RHS of the same function and I do not know how to do that.

In terms of the Euler's method, it is a explicit method so the formula is given by:

$Y_{n+1}=Y_n+hf(t_n,Y_n)$ which is already in the explicit formula and for my function I will obtain:

$Y_{n+1}=Y_n+h(\frac{8 \cos(4t_n)+y_n}{4})$

For the 2nd order Taylor method I really do not understand this method.

Can anyone help me on this and point me to a website where I can find more information about these methods?

Thanks

Best Answer

Both Explicit Trapezoid method and Heun's method are 2nd order Runge Kutta methods in this form:

$$y_{n+1}=y_n+b_1hf(t_n,y_n)+b_2hf(t_i+c_2h,y_n+a_{21}hf(t_n,y_n)).$$

For them to be second order, they need to satisfy: $$b_1+b_2=1, b_2c_1=1/2, a_{21}b_2=1/2.$$

For Explicit Trapezoid, we have

$$b_1=b_2=1/2, c_2=1, a_{21}=1.$$

For Heun's, we have $$b_1=1/4, b_2=3/4, c_2=2/3, a_{21}=1.$$

For your example, Heun's method gives:

$$y_n+\frac{h}{4}\cdot\frac{8\cos(4t_n)+y_n}{4}+\frac{3h}{4}\frac{8\cos(4(t_n+\frac{2}{3}h))+y_n+h\cdot\frac{8\cos(4t_n)+y_n}{4}}{4}.$$

I think you can find this in any Numerical Analysis book.

Your Trapezoid method formula is correct, but you didn't apply it correctly in your example. It should be:

$$y_n+\frac{h}{2}\frac{8\cos(4t_n)+y_n}{4}+\frac{h}{2}\frac{8\cos(4(t_n+h))+h\frac{8\cos(4t_n)+y_n}{4}}{4}.$$

The 2nd order Taylor method is a totally different method. It uses Taylor expansion, so requires more partial derivatives of $f$:

$$y_{n+1}=y_n+hf(t_n,y_n)+\frac{h^2}{2}\left(\frac{\partial f}{\partial t}(t_n,y_n)+\frac{\partial f}{\partial y}(t_n,y_n)f(t_n,y_n)\right).$$

Applying to your example, we get $$y_{n+1}=y_n+h\frac{8\cos(4t_n)+y_n}{4}+\frac{h^2}{2}\left(-8\sin(4t_n)+\frac{1}{4}\frac{8\cos(4t_n)+y_n}{4}\right).$$

Related Solutions

[Math] Rewrite matrix equation for Euler method and Improved Euler method

Euler's Method works as follows: Given an ODE $y'(t) = f(t,y(t)$ we can compute $y_{n+1} = y_n +h f(t_n,y(t_n)) $ with some given stepsize $h$. In your case ( even though you are speaking of $y_{n+1}$, but your ODE reads $x'=Ax+g$, I assume you then have the ODE $y'=Ay+g$ ) $y_{n+1} = y_n +h (A y_n +g) = (A+h)y_n+ hg $.

Improved Euler method can be split into two equation, but it should rather be called two steps. You basically perform one Euler step, but this time you call the solution $y^*_{n+1}$. You then use the average of the slope at $y_n$ and $y^*_{n+1}$ to approximate the slope for your step $y_{n+1} = y_{n} +h\cdot s$ with $s = (f(y_n)+ f(y_{n+1}^*)/2$. Applied to your problem this would result in the following term. \begin{align}y_{n+1} &=y_n + h\cdot(f(y_n)+ f(y_{n+1}^*)/2 \\ &= y_n+ \frac{h}{2}( A y_n+g+(A+h)y_n+hg) \text{ by Euler method for }y^*_{n+1} \\ &= (1+{Ah}+ \frac{ h^2}{2})y_n +(\frac{h}{2}+\frac{h^2}{2})g\end{align}

Thomas

[Math] Stability and convergence for Heun’s method

Stability of explicit Runge-Kutta methods around a fixed point or a slow moving equilibrium is obtained when $Lh\sim 1$ with $L$ the $y$-Lipschitz constant and $h$ the step size. The exact constant on the right side depends on the order of the method, usually it is between 1 and 5. For first and second order methods that is only a guideline as the stability region does not contain a semi-circle around the origin in the negative complex half-plane.

To get a practical feeling for this, just solve your problem for various step sizes over the integration interval from $0$ to $1/3$, as $y(1/3)=0.1e^{-50}$ is a really small number, relative to the initial value and floating point precision.

def Heun(f,t0,y0,tf,h):
    def HeunStep(t,y,h): 
        k1 = h*f(t,y); 
        k2 = h*f(t+h,y+k1); 
        return y+0.5*(k1+k2);
    t, y = t0, y0
    while t+1.1*h<tf: t, y = t+h, HeunStep(t,y,h)
    return HeunStep(t,y,tf-t)

def f(t,y): return 49-150*(t+y)
t0,y0 = 0.0,1./3+0.1
for h in np.linspace(3.5,0.5,16)/150:
    print "h=%.4e 150*h=%.4e, y(1/3)=%.12g"%(h,150*h, Heun(f,t0,y0,1./3,h))

with the results

h=2.3333e-02 150*h=3.5000e+00, y(1/3)=3382845.45109
h=2.2000e-02 150*h=3.3000e+00, y(1/3)=1820496.14682
h=2.0667e-02 150*h=3.1000e+00, y(1/3)=558708.429462
h=1.9333e-02 150*h=2.9000e+00, y(1/3)=79763.1519923
h=1.8000e-02 150*h=2.7000e+00, y(1/3)=9204.09082254
h=1.6667e-02 150*h=2.5000e+00, y(1/3)=1648.41784102
h=1.5333e-02 150*h=2.3000e+00, y(1/3)=37.6118108918
h=1.4000e-02 150*h=2.1000e+00, y(1/3)=0.740437945724
h=1.2667e-02 150*h=1.9000e+00, y(1/3)=0.00432803787737
h=1.1333e-02 150*h=1.7000e+00, y(1/3)=1.0688786734e-05
h=1.0000e-02 150*h=1.5000e+00, y(1/3)=1.14794370277e-08
h=8.6667e-03 150*h=1.3000e+00, y(1/3)=5.57823318786e-12
h=7.3333e-03 150*h=1.1000e+00, y(1/3)=2.8015784137e-15
h=6.0000e-03 150*h=9.0000e-01, y(1/3)=2.34187669257e-17
h=4.6667e-03 150*h=7.0000e-01, y(1/3)=-1.38777878078e-17
h=3.3333e-03 150*h=5.0000e-01, y(1/3)=2.60208521397e-17

As the stability region suggests, $0<h<\frac2{150}$ are valid choices in the table above, giving results close to zero. Sensible results that are practially zero as floating point numbers (and relative to the initial value) are only obtained for about $h<\frac{1.5}{150}$.

Best Answer

Related Solutions

[Math] Rewrite matrix equation for Euler method and Improved Euler method

[Math] Stability and convergence for Heun’s method

Related Question