[Numerical] Euler Backward method

[Nick89]

Hi,

I need to solve the following simple ODE with both the Euler Forward and Euler Backward numerical methods. I also need to answer for which values of T this can still be calculated:
$$y'(t) = \frac{-1}{2y(t)}$$
$$y(0) = 2$$
$$t \in (0, T]$$

Obviously the analytical solution is
$$y(t) = \sqrt{4 - t}$$
So it would seem T must be between 0 and 4 for the root to be real.The Euler Forward method, if you are not familiar with it, follows from doing a taylor expansion and reads:
$$y_{n+1} = y_n + h y'_n$$
Here, h is the 'step size'.
This method is called the explicit method because y_n+1 can be calculated explicitly by simply evaluating the right hand side, which is known.
If an initial value is known, any next value can be calculated from this. Because y(0) is given (=2) I have been able to put this into a Matlab script and graph the numerical solution, which matched the analytical solution perfectly as far as I could tell.Next, I also need to graph the solution using the Euler Backward method, which looks similar but is not quite the same:
$$y_n = y_{n+1} - h y'_{n+1}$$
This method is also called the implicit method because y_n+1 cannot be calculated explicitly by evaluating the right hand side. Instead it needs to be solved for y_n+1.

Unfortunately, I have no clue how to build 'solve for y_n+1' into Matlab, and I'm pretty sure we are not required to do this.
So I tried just entering the ODE to get:
$$y_n = y_{n+1} + \frac{h}{2 y_{n+1}}$$

Solving this for y_n+1 by hand, I obviously get a quadratic equation for y_n+1, for which there are two solutions:
$$y_{n+1} = \frac{ -y_n \pm \sqrt{y_n^2 - 2h}}{2}$$Now finally my question... Of course I can use the equation above to evaluate y_n+1, but how do I know whether to use the + or the - ?
I really have no idea..? I could use the + all the time, or the - all the time, but I have no reason to believe that I should not use +, -, +, - etc or something like that. It does not sound logical to require it to use +,-,+,- etc but still, I simply don't know...

How can I determine whether to use the + or -?

Thanks!

 

[matematikawan]

Just wondering. What advantages does the Euler Backward method has over the Euler Forward?

 

[D H]

Solving this for y_n+1 by hand, I obviously get a quadratic equation for y_n+1, for which there are two solutions:
$$y_{n+1} = \frac{ -y_n \pm \sqrt{y_n^2 - 2h}}{2}$$

Note well: You have a sign error. The two solutions are

$$y_{n+1} = \frac{y_n \pm \sqrt{y_n^2 - 2h}}{2}$$

Which to choose? First, it helps to rewrite the above as

$$y_{n+1} = y_n\left(\frac{1 \pm \sqrt{1 - 2h/y_n^2}}{2}\right)$$

Examine what happens as the step size becomes very small. One choice makes a lot of sense. The other choice is complete nonsense.

Just wondering. What advantages does the Euler Backward method has over the Euler Forward?

Note well: Euler techniques almost always yield very poor results.

That said, the advantage of using implicit integration techniques is stability (but typically at the cost of increased complexity and sometimes decreased accuracy). A simple example is exponential decay with a time constant $\tau$:

$$\frac{dx}{dt} = -{x(t)}/{\tau}$$

Here $\tau>0$. Forward Euler integration with step size $\Delta t$ leads to

$$x_{n+1} = x_n(1 - {\Delta t}/{\tau})$$

If you choose a time step that is larger than the time constant the above will diverge. Even if the time step is less than the time constant, the results can be quite bad. You choose a time step that is much less than the time constant.

In comparison, backwards Euler leads to

$$x_{n+1} = \frac{x_n}{1 + {\Delta t}/{\tau}}$$

which is always stable, even with a bad choice of a step size.

 

[Nick89]

Note well: You have a sign error. The two solutions are

$$y_{n+1} = \frac{y_n \pm \sqrt{y_n^2 - 2h}}{2}$$

Oops, I think you're right.

Which to choose? First, it helps to rewrite the above as

$$y_{n+1} = y_n\left(\frac{1 \pm \sqrt{1 - 2h/y_n^2}}{2}\right)$$

Examine what happens as the step size becomes very small. One choice makes a lot of sense. The other choice is complete nonsense.

When I use the + sign it converges to y_n, when I use the - sign it becomes zero. That makes sense, thanks.

I ended up using the Newton iteration method to solve the equation numerically (just because it is more general and can be applied to other functions as well), but it yields the same result (just a little slower).

Thanks for the help!

 

[matematikawan]

Note well: Euler techniques almost always yield very poor results.

That said, the advantage of using implicit integration techniques is stability (but typically at the cost of increased complexity and sometimes decreased accuracy). A simple example is exponential decay with a time constant $\tau$:

$$\frac{dx}{dt} = -{x(t)}/{\tau}$$

Here $\tau>0$. Forward Euler integration with step size $\Delta t$ leads to

$$x_{n+1} = x_n(1 - {\Delta t}/{\tau})$$

If you choose a time step that is larger than the time constant the above will diverge. Even if the time step is less than the time constant, the results can be quite bad. You choose a time step that is much less than the time constant.

In comparison, backwards Euler leads to

$$x_{n+1} = \frac{x_n}{1 + {\Delta t}/{\tau}}$$

which is always stable, even with a bad choice of a step size.

This interesting. I have only basic knowledge of Euler's method. Never been introduced to the term forward or backward method before.

Of course in practise people seldom use Euler forward method. It is only for illustration of numerical DE mthod. I do not know about Euler backward method.

So it is all about stability in using Euler backward method. But it look like that this method only work if we could algebraically solve y_n=y_n+1-hy'_n+1 for y_n+1. That a big problem. Not many DE can be solve with this method.


Another curious question. Do they also have such thing as Modified Euler Backward Method?

 

[Nick89]

But it look like that this method only work if we could algebraically solve y_n=y_n+1-hy'_n+1 for y_n+1. That a big problem. Not many DE can be solve with this method.

Actually, you can just enter the DE for y'_n+1 which will probably be a function of y_n+1. The 'problem' is now reduced simply to solving a (possibly non-linear) equation. Many equations like this can be solved analytically, and if they cannot there are always numerical methods.

 

[matematikawan]

Still need to be convinced on using Euler backward method. Please explain how to solve this IVP.

y' = sin y subject to let say y(0)=1.

We start with this equation
y_n = y_n+1 - h sin(y_n+1)

Then what next? It is not easy to express y_n+1 in term of y_n.

 

[D H]

The typical approach to using implicit methods for non-linear ODEs is to linearize the derivative function, y'=f(x,y), where x is the independent variable. The method is no longer purely implicit and is better called semi-implicit. Implicit Euler integration uses

$$y_{n+1}=y_n+h f(x_{n+1},y_{n+1})$$

Suppose the derivative function does not depend on the dependent variable. Linearizing the derivative function about y_n,

$$f(y_{n+1}) \approx f(y_n) + (y_{n+1}-y_n)\left.\frac{\partial f}{\partial y}\right|_{y_n}$$

With this, the semi-implicit Euler integration scheme becomes

$$y_{n+1} = y_n + h \left( f(y_n)+(y_{n+1}-y_n)\left.\frac{\partial f}{\partial y}\right|_{y_n} \right)$$

Solving for y_n+1,

$$y_{n+1} = y_n + h\left(1-h\left.\frac{\partial f}{\partial y}\right|_{y_n}\right)^{-1} f(y_n)$$

In the problem at hand, y'=sin y, this becomes

$$y_{n+1} = y_n + h\left(1-h\cos y_n\right)^{-1} \sin y_n$$

 

[matematikawan]

Thanks DH for the explanation.

So there can be two type of Euler Backward Method.
(i) the implicit y_n=y_n+1 - h y'_n+1 and
(ii) the semi-implicit $$y_{n+1}=y_n+h f(x_{n+1},y_{n+1})$$ .

Or is there others?

I hope they have the theorem which state that this method is always stable. I will try later to write a MATLAB script later to see what result this method gives. Has anyone seen the analytic solution for
$$\frac{dy}{dx} = \sin(y)$$
before?

 

[matematikawan]

I have written the MATLAB script. If only I know the analytic solution to compare.


% Solving IVP y'=sin(y) , y(0)=1
% Using Euler Backward Method

h=0.1; % step size
y(1)=1; % MATLAB index start from 1. y0=y(1)
xi=0; xf=10; % domain [xi , xf] num_pt=(xf-xi)/h+1 = 101
for n=2:101
y(n)=y(n-1)+h*sin(y(n-1))/(1-h*cos(y(n-1)));
end
x=xi:h:xf;
plot(x,y)

 

[matematikawan]

OK got it.


$$\int \csc (y) dy = x + c$$
$$\log (\tan (\frac{y}{2})) = x + c$$

ref: http://integrals.wolfram.com/index.jsp?expr=Csc[x]&random=false

$$y = 2\tan^{-1} (\tan (0.5) e^x)$$

 

[ametista82]

I have the same problem
I have this force...
F= k*l-ln/ln*ppj/||p*pj|| where p and pj is the position of particle p and pj, l and ln are the length of a spring. how can use the implicit integration in this ode?
and If i have a costant force that depends only by an constant acceleration? thank you for your help

 

[matematikawan]

Not sure whether I could still remember this stuff. If you could please write down the related differential equation (if possible in latex) and initial condition it will be more clearer your problem.

 

[ametista82]

First: thank for your answer!
Sorry...I read a Baraff's paper about implicit integration..(sorry but I don't use Latex..)...so by this paper:

deltax = h(v0 +deltav) where deltax=x(t0+h)-x(t0)

(I -hdf/dv -h^2df/dx)deltav=h(f0+hdf/dxv0)

so the problem is to solve for deltav, so if my force is

F=fg + k* l-ln/ln ppj/|ppj| where fg is gravitational force, l is the natural length of the spring and ppj are the position of two particle betwen the string.

So how can use implicit integration? And another question if I have only a costant force as F=ma how can use an implicit integration?

 

[klunkerz]

Solving this for y_n+1 by hand, I obviously get a quadratic equation for y_n+1, for which there are two solutions:
$$y_{n+1} = \frac{ -y_n \pm \sqrt{y_n^2 - 2h}}{2}$$



Hi,
How to solve the above equation by hand?
Can anyone please explain it to me?

 

[matematikawan]

I though everybody has access to at least a calculator to compute square root of a number.

Referring back to DH post (the 3rd post)
$$y_{n+1} = y_n\left(\frac{1 + \sqrt{1 - 2h/y_n^2}}{2}\right)$$

Apply the binomial expansion
$$y_{n+1} = \frac{y_n}{2}\left(1 + 1-\frac{2h}{2y_n^2}+...\right)$$
$$\approx y_n\left(1 -\frac{h}{2y_n^2}\right)$$

which can be calculated using a 'slide rule'/table