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

Backward Euler / Implicit Integration Implementation

  1. Nov 13, 2009 #1
    Hi all,

    I am trying to implement the backward euler integration (in c++) for the pendulum problem. I have the forward euler implemented, but frankly I don't know where to even start from for the implicit integration. I understand the update expressions for implicit, and of course the pendulum motion (theta'' = -gravity / length * sin (theta) ).

    Almost all content online explaining implicit integration give the update expressions and say that you have to solve for x (t+1), but I haven't found any resources that actually explain how to solve them.

    I am not looking for libraries that do this already, I would like to learn the process and do it myself (except maybe for matrix inversion etc). If you understand the process of solving for x(t+1) would you please enlighten me? Do I need to find the derivative of the acceleration? Do I need to write the motion function ( -gravity / length * sin (theta) ) in matrix notation?

    Very confused.

    Thank you so much!
  2. jcsd
  3. Nov 14, 2009 #2
    For forward Euler I had seen a number of times similar problem posted to this forum. Normally people use RK4 method to solve.

    Just let [tex] x_1=\theta , x_2=\dot{\theta}[/tex]
    then expressed your equation as a system of DE
    [tex]\dot{\vec x} = \boldsymbol A \vec x[/tex]

    where A=[0 1; -g/L*sin(x1) 0]
  4. Nov 14, 2009 #3

    Thanks, but that's for forward euler? I'm asking about backward euler.
  5. Nov 14, 2009 #4
    Yeap I didn't answer your question and I know that. I'm trying to express what little knowledge that I have on the subject :biggrin:

    Now I remember that I had seen a similar problem, but much simplier in the forum.

    Looking back at that thread I think, for the Backward Euler we have to work on semi-implicit (I didn't create this term). My guess for the iteration is

    Xn+1 = Xn + h*(I - hJ)-1*F(Xn)

    where X=[x1 x2]t
    J is the Jacobian matrix [0 1; -g/L*cos(x1) 0]
    F(X)= [x2 -g/L*sin(x1)]t
  6. Dec 14, 2009 #5
    Hi! I have the same problem! Ihave to implement some natural effects with implicit integration,for example if I have a particle with an acceleration toward a point (orbit) so there is this force:

    F(m, epsilon)=m/(r^2)+epsilon
    with euler explicit I have this code:

    Code (Text):
    loat magdt = magnitude * dt;
    float max_radiusSqr = max_radius * max_radius;
    if(max_radiusSqr < P_MAXFLOAT) {
    for (ParticleList::iterator it = ibegin; it != iend; it++) {
    Particle_t &m = (*it);
    pVec dir(center - m.pos// Step velocity with acceleration
    if(rSqr < max_radiusSqr)
    m.vel += dir * (magdt / (sqrtf(rSqr) * (rSqr + epsilon)));//sqrtf is for to normalize rSqr

    How can I implement this with implicit integration???

    thank you..
  7. Dec 16, 2009 #6
    I'm not that familier with C++ code. If you could express your problem mathematically, probably we could come out with the algorithm.
  8. Dec 16, 2009 #7
    ok thank you! sorry!
    My problem is:
    I have a particle systems and I have to simulate an effect called atom, and with Explicit Euler for to calculate particle's acceleration, I have this:

    new_vel = old_vel + dt * (m /( r^2+epsilon));

    epsilon serves to dampen the inverse square so that f(r) does not approach infinity as r approches zero, r is the distance from gravity source.
    Now my question is:how can I calcute the same new_vel using an implicit integration?
  9. Dec 17, 2009 #8
    I'm not sure whether I properly understand your problem here. I think your r is also a dependant variable in t. Epsilon and m are constants.

    Vn+1 = Vn + dt*f(rn+1,tn+1)

    We need another DE in dr/dt.
Share this great discussion with others via Reddit, Google+, Twitter, or Facebook