Implicit Euler Integration - Size of Stability Region

Calculus Level pending

Consider the formula for implicit Euler integration of a complex exponential. Notice that the only difference between this expression and the expression from the earlier problem is that the derivative term is from the present rather than the past.

y ( t ) = y ( t Δ t ) + y ˙ ( t ) Δ t y ( t ) = e λ t y ˙ ( t ) = λ y ( t ) λ = a + j b j = 1 y(t) = y(t - \Delta t) + \dot{y}(t) \, \Delta t \\ y(t) = e^{\lambda \, t} \\ \dot{y}(t) = \lambda \, y(t) \\ \lambda = a + j \, b \\ j = \sqrt{-1}

The complex exponential is processed sequentially in time with a discrete time interval Δ t \Delta t . In order for the algorithm to be stable (non-divergent) for a given Δ t \Delta t , the following condition must hold true:

y ( t ) y ( t Δ t ) 1 \large{\Big| \frac{y(t)}{y(t - \Delta t)} \Big | \leq 1}

In other words, the modulus of the ratio of the present value of y y to the previous value of y y must be less than or equal to 1 1 . Note that this ratio will generally be a complex number.

Let's examine the case where Δ t = 1 \Delta t = 1 . Suppose we were to plot all non-divergent ( a , b ) (a,b) points on a two-dimensional plane, with the horizontal axis corresponding to the a a values and the vertical axis corresponding to the b b values. This is known as the region of stability.

Is the area of the stability region finite or infinite? In other words, is there a finite number larger than the area of the stability region?

Infinite Finite

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

Chris Lewis
May 27, 2019

Substituting y ˙ ( t ) = λ y ( t ) \dot{y}(t)=\lambda y(t) into the first equation gives

y ( t ) = y ( t Δ t ) + λ y ( t ) Δ t y(t)=y(t-\Delta t)+\lambda y(t) \Delta t

Rearranging,

y ( t ) y ( t Δ t ) = 1 1 λ Δ t \frac{y(t)}{y(t-\Delta t)}=\frac{1}{1-\lambda \Delta t}

Since λ = a + b i \lambda = a+bi and Δ t = 1 \Delta t=1 , we require the region of C \mathbb{C} satisfying

1 a 1 + b i 1 \frac{1}{\left|a-1+bi \right|} \le 1

or

a 1 + b i 1 \left|a-1+bi \right| \ge 1

Clearly this region is infinite; specifically, it is everywhere in the complex plane except the interior of the circle centred at 1 1 with radius 1 1 .

Thanks for the solution. Actually, the instability region is centered at + 1 +1 with radius 1 1 . So the stability region contains the entire left-hand plane, meaning that the algorithm converges for any decaying exponential, and even converges for many exponentially growing inputs.

Steven Chase - 2 years ago

Log in to reply

Yes, of course that's the centre, silly me. Typo corrected, thank you!

Really enjoying these illuminating problems, by the way, especially how you intersperse these more theoretical ones with actual applications. It's a while since I studied these so good to get a refresher course.

One query: is the method stable everywhere on the circle itself? It seems bounded, but infinite, orbits would form, and not necessarily converge. Though I may just be getting muddled here.

Chris Lewis - 2 years ago

Log in to reply

Thanks, glad you are liking them. Exactly on the boundary, the following holds true:

y ( t ) y ( t Δ t ) = e j θ \large{\frac{y(t)}{y(t - \Delta t)} = e^{j \theta}}

So in this case, y ( t ) y(t) is either a stationary complex number (if θ = 0 \theta = 0 ), or a vector of constant length rotating in the complex plane (if θ 0 \theta \neq 0 ). In the second scenario, the projection of the rotating vector onto one axis yields a sinusoid, which connects us back to the Euler equation.

But yes, you're right. On the boundary, we have neither growth nor decay.

Steven Chase - 2 years ago

0 pending reports

×

Problem Loading...

Note Loading...

Set Loading...