Dismiss Notice
Join Physics Forums Today!
The friendliest, high quality science and math community on the planet! Everyone who loves science is here!

[Numerical] Euler Backward method

  1. Mar 23, 2009 #1
    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:
    [tex]y'(t) = \frac{-1}{2y(t)}[/tex]
    [tex]y(0) = 2[/tex]
    [tex]t \in (0, T][/tex]

    Obviously the analytical solution is
    [tex]y(t) = \sqrt{4 - t}[/tex]
    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:
    [tex]y_{n+1} = y_n + h y'_n[/tex]
    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:
    [tex]y_n = y_{n+1} - h y'_{n+1}[/tex]
    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:
    [tex]y_n = y_{n+1} + \frac{h}{2 y_{n+1}}[/tex]

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


    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!
     
  2. jcsd
  3. Mar 23, 2009 #2
    Just wondering. What advantages does the Euler Backward method has over the Euler Forward?
     
  4. Mar 24, 2009 #3

    D H

    User Avatar
    Staff Emeritus
    Science Advisor

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

    [tex]y_{n+1} = \frac{y_n \pm \sqrt{y_n^2 - 2h}}{2}[/tex]

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

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

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

    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 [itex]\tau[/itex]:

    [tex]\frac{dx}{dt} = -{x(t)}/{\tau}[/tex]

    Here [itex]\tau>0[/itex]. Forward Euler integration with step size [itex]\Delta t[/itex] leads to

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

    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

    [tex]x_{n+1} = \frac{x_n}{1 + {\Delta t}/{\tau}}[/tex]

    which is always stable, even with a bad choice of a step size.
     
  5. Mar 24, 2009 #4
    Oops, I think you're right.


    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!
     
  6. Mar 24, 2009 #5
    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 yn=yn+1-hy'n+1 for yn+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?
     
  7. Mar 25, 2009 #6
    Actually, you can just enter the DE for y'n+1 which will probably be a function of yn+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.
     
  8. Mar 26, 2009 #7
    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
    yn = yn+1 - h sin(yn+1)

    Then what next? It is not easy to express yn+1 in term of yn.
     
  9. Mar 26, 2009 #8

    D H

    User Avatar
    Staff Emeritus
    Science Advisor

    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

    [tex]y_{n+1}=y_n+h f(x_{n+1},y_{n+1})[/tex]

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

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

    With this, the semi-implicit Euler integration scheme becomes

    [tex]
    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)
    [/tex]

    Solving for yn+1,

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

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

    [tex]
    y_{n+1} = y_n + h\left(1-h\cos y_n\right)^{-1} \sin y_n
    [/tex]
     
  10. Mar 26, 2009 #9
    Thanks DH for the explaination.

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

    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
    [tex]\frac{dy}{dx} = \sin(y)[/tex]
    before?
     
  11. Mar 29, 2009 #10
    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)
     

    Attached Files:

    • ebm.jpg
      ebm.jpg
      File size:
      8.5 KB
      Views:
      225
  12. Mar 29, 2009 #11
  13. Feb 4, 2010 #12
    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
     
  14. Feb 4, 2010 #13
    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.
     
  15. Feb 4, 2010 #14
    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?
     
  16. Nov 5, 2010 #15
    Solving this for y_n+1 by hand, I obviously get a quadratic equation for y_n+1, for which there are two solutions:
    [tex]y_{n+1} = \frac{ -y_n \pm \sqrt{y_n^2 - 2h}}{2}[/tex]



    Hi,
    How to solve the above equation by hand?
    Can anyone please explain it to me?
     
  17. Nov 6, 2010 #16
    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)
    [tex] y_{n+1} = y_n\left(\frac{1 + \sqrt{1 - 2h/y_n^2}}{2}\right) [/tex]

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

    which can be calculated using a 'slide rule'/table :wink:
     
Know someone interested in this topic? Share this thread via Reddit, Google+, Twitter, or Facebook