Numerical Stability of Euler Variants

Calculus Level 3

Consider the classic explicit Euler numerical integration technique:

x k = x k − 1 + x ˙ k − 1 Δ t \large{x_k = x_{k-1} + \dot{x}_{k-1} \, \Delta t}

And a modified version:

x k = x k − 1 + x ˙ k − 1 Δ t + x ¨ k − 1 ( Δ t ) 2 \large{x_k = x_{k-1} + \dot{x}_{k-1} \, \Delta t + \ddot{x}_{k-1} \, (\Delta t)^2 }

Suppose we use these to discretely approximate the following continuous-time signal:

x ( t ) = e α t \large{x(t) = e^{\alpha t}}

In the above equation, α \alpha is a complex number. Define "stability" as the following being true for all processing steps:

∣ x k x k − 1 ∣ < 1 \large{\Big| \frac{x_k}{x_{k-1}} \Big| < 1}

For each integration method, determine the region of stability in the complex α \alpha plane, as well as the area of each region.

For Δ t = 1 \Delta t = 1 , what is the following ratio?:

Area of modified Euler stability region Area of explicit Euler stability region \frac{\text{Area of modified Euler stability region}}{\text{Area of explicit Euler stability region}}

Inspiration

Bonus: Make scatter plots of both stability regions. Actually, this is the best way to map out the stability regions anyway


The answer is 0.8394.

Steven Chase by Steven Chase , 37, USA

This section requires Javascript.
You are seeing this because something didn't load right. We suggest you, (a) try refreshing the page, (b) enabling javascript if it is disabled on your browser and, finally, (c) loading the non-javascript version of this page . We're sorry about the hassle.

1 solution

Mark Hennings
Nov 7, 2019

Since x ˙ = α x \dot{x} = \alpha x and x ¨ = α 2 x \ddot{x} = \alpha^2x , the explicit Euler formula in this case is simply x k = ( 1 + α ) x k − 1 x_k \; = \; (1 + \alpha)x_{k-1} and hence the stability condition is ∣ 1 + α ∣ < 1 |1+\alpha| < 1 . Thus α \alpha lies in the disk with centre − 1 -1 and radius 1 1 , and hence the explicit stability region has area π \pi .

On the other hand, the modified formula is simply x k = ( 1 + α + α 2 ) x k − 1 x_k \; = \; (1 + \alpha + \alpha^2)x_{k-1} which means that the stability condition is ∣ ( α + 1 2 ) 2 + 3 4 ∣ = ∣ α 2 + α + 1 ∣ < 1 \left|\big(\alpha+\tfrac12\big)^2 + \tfrac34\right| = |\alpha^2 + \alpha + 1| < 1 which means that ( α + 1 2 ) 2 (\alpha+\tfrac12)^2 must lie in the circle with centre − 3 4 -\tfrac34 and radius 1 1 . In polar coordinates this means that ( α + 1 2 ) 2 = r e i θ (\alpha + \tfrac12)^2 = re^{i\theta} where 0 ≤ r < 1 4 ( 7 + 9 cos ⁡ 2 θ − 3 cos ⁡ θ ) − π ≤ θ ≤ π 0 \le r < \tfrac14\left(\sqrt{7 + 9\cos^2\theta} - 3\cos\theta\right) \hspace{2cm} -\pi \le \theta \le \pi Thus we deduce that α = r e i θ − 1 2 \alpha = re^{i\theta} - \tfrac12 where 0 ≤ r < R ( θ ) = 1 2 7 + 9 cos ⁡ 2 2 θ − 3 cos ⁡ 2 θ − π ≤ θ ≤ π 0 \le r < R(\theta) = \tfrac12\sqrt{\sqrt{7 + 9\cos^22\theta} - 3\cos2\theta} \hspace{2cm} -\pi \le \theta \le \pi A graph of the modified stability region is shown. Its area is 4 ∫ 0 1 2 π 1 2 R ( θ ) 2 d θ = 1 2 ∫ 0 1 2 π ( 7 + 9 cos ⁡ 2 2 θ − 3 cos ⁡ 2 θ ) d θ = 2 E ( 9 16 ) 4\int_0^{\frac12\pi} \tfrac12R(\theta)^2\,d\theta \; = \; \tfrac12\int_0^{\frac12\pi}\left(\sqrt{7 + 9\cos^22\theta} - 3\cos2\theta\right)\,d\theta \; = \; 2E(\tfrac{9}{16}) where E E is the complete elliptic integral. This makes the ratio of areas 2 π E ( 9 16 ) = 0.839365 \tfrac{2}{\pi}E(\tfrac{9}{16}) = \boxed{0.839365} . No scatter diagrams required.

Is the term "stability" suitable in this second case? As can be seen, the "stability region" contains numbers with positive real part, and hence for which ∣ x ( k ) ∣ |x(k)| will diverge to ∞ \infty as k → ∞ k \to \infty . Since we are trying to approximate x ( k ) x(k) by x k x_k , we don't want ∣ x k + 1 ∣ < ∣ x k ∣ |x_{k+1}| < |x_k| in this case, do we?

2 Helpful 2 Interesting 2 Brilliant 0 Confused

Indeed, by "stable" I just mean "doesn't blow up". As you say, there are some exponentially increasing signals that will exponentially decrease when numerically approximated. So the results for those cases are bogus. Similarly, there are exponentially decaying signals that will exponentially grow when numerically approximated. The trapezoidal method, by contrast, contains the entire left plane and nothing in the right plane.

Steven Chase - 1 year, 7 months ago

Log in to reply

True. The bottom line is that all numerical techniques are (eventually) inaccurate, and that errors accumulate. The question is how fast the methods diverge from the actual value, enabling a cost-benefit analysis conversation about how small h = Δ t h = \Delta t must be to enable an acceptable level of error, given that a smaller value of h h will require more iterations. The trapezoidal technique has a cubic local truncation error, so the rate of divergence from the actual value is proportional to h 3 h^3 instead of h h ; this makes it better than other schemes, but still inaccurate if pushed too far.

Mark Hennings - 1 year, 7 months ago

0 pending reports

×

Problem Loading...

Note Loading...

Set Loading...