Numerical integration - Modified Euler

Calculus Level 3

Consider the differential equation of the following linear oscillator:

$\boxed{\ddot{x} +100x = 0} \ ; \ x(0) = 1 \ ; \ \dot{x} = 0$

Compute the analytical solution to this equation. Following this, numerically integrate the above equation using the following two methods:

Method 1: Explicit Euler $\dot{x}(k+1) =\dot{x}(k) + h\ddot{x}(k)$ $x(k+1) =x(k) + h\dot{x}(k)$

Method 2: Modified Euler $\dot{x}(k+1) =\dot{x}(k) + h\ddot{x}(k)$ $x(k+1) =x(k) + h\dot{x}(k+1)$

Here $k$ denotes the time index. Numerical integrate for two seconds. Let the analytical solution be $x_a$ , the explicit Euler solution be $x_{exp}$ and the modified Euler solution be $x_{mod}$ . At the time instant $t = 2$ compute the following expressions:

$\mid x_{mod} - x_a \mid$ $\mid x_{exp} - x_a \mid$

Enter the result with a smaller magnitude as your answer. Take the time step $h = 0.01$ .

Bonus: Which numerical method would you prefer and why? Also, inspect the derivative signals of $x$ and state your conclusions.

The answer is 0.053504.

1 solution

Steven Chase
Nov 6, 2019

Simulation code and results are below. The modified Euler technique is more accurate because the position is updated based on the first and second derivatives from the previous time step, whereas the explicit Euler position is updated only with the first derivative. For the given time step, the performance difference is dramatic (see plot).

Note that for small time steps (ex: $10^{-5}$ ), all three methods agree closely. In practice, almost everything I program is for non-real-time applications in which speed is not important. So I generally prefer to run explicit Euler with a very small time step, due to the simplicity of the method. To check the results, I typically do multiple runs with time steps spanning two or three orders of magnitude, and if the results are the same each time, I am satisfied with the answer. In real-time or performance-sensitive applications, we could consider more sophisticated integration methods.

Another thing I want to look at is the numerical stability of the modified Euler method given the input signal $x = A e^{ \alpha t}$ . Given that input signal and a particular time step, what is the $\alpha$ region in the complex plane for which the method is stable? And how does that region compare to the explicit Euler stability region?

import math

tf = 2.0
dt = 0.01

#####################################

t = 0.0

xe = 1.0 
xed = 0.0
xedd = -100.0*xe

xm = 1.0
xmd = 0.0
xmdd = -100.0*xm

xa = math.cos(10.0*t)

print "t xa xe xm xed xmd xedd xmdd"

while t <= tf:

    print t,xa,xe,xm,xed,xmd,xedd,xmdd

    t = t + dt

    xa = math.cos(10.0*t)    # Analytical value

    xe = xe + xed*dt         # Explicit Euler
    xed = xed + xedd * dt
    xedd = -100.0*xe

    xm = xm + xmd*dt + xmdd*(dt**2.0)      # Modified Euler
    xmd = xmd + xmdd * dt
    xmdd = -100.0*xm

#####################################

print ""
print ""

print xa
print xe
print xm
print ""
E1 = math.fabs(xm - xa)
E2 = math.fabs(xe - xa)

print E1
print E2

#0.408082061813
#1.26488581312
#0.354578224809

#0.0535038370045
#0.856803751308

@Karan Chatrath Mind if I post the stability follow-up to this?

Steven Chase - 1 year, 7 months ago

Thanks for the solution. I don't mind. Looking forward to the follow up

Karan Chatrath - 1 year, 7 months ago

It's up now. By the way, which would you choose as the favorite?

Steven Chase - 1 year, 7 months ago

@Steven Chase – Thanks. I will attempt it later today.

While I do appreciate the simplicity of explicit Euler, I prefer the modified algorithm. It is higher-order accurate and it is not very hard to program as well. The plot in the solution shows a remarkable difference between the two algorithms. This observation is primarily why I chose to make a problem out of it.

Karan Chatrath - 1 year, 7 months ago

For oscillating systems such as this, I would suggest the Runge-Kutta family of methods are amazing. The problem with the Explicit Euler method is extremely apparent in oscillatory systems, where the energy increases very visibly, even with smaller time-steps.

Krishna Karthik - 1 year, 1 month ago

In that way, every single numerical solving method for differential equations is better than Explicit Euler!

Krishna Karthik - 1 year, 1 month ago

Numerical integration - Modified Euler

The answer is 0.053504.

1 solution

0 pending reports