Backward Euler / Implicit Integration Implementation

[jacki]

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!

 

[matematikawan]

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 x_1=\theta , x_2=\dot{\theta}
then expressed your equation as a system of DE
\dot{\vec x} = \boldsymbol A \vec x

where A=[0 1; -g/L*sin(x1) 0]

 

[jacki]

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 x_1=\theta , x_2=\dot{\theta}
then expressed your equation as a system of DE
\dot{\vec x} = \boldsymbol A \vec x

where A=[0 1; -g/L*sin(x1) 0]

Thanks, but that's for forward euler? I'm asking about backward euler.

 

[matematikawan]

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

Now I remember that I had seen a similar problem, but much simplier in the forum.
https://www.physicsforums.com/showthread.php?t=301890

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

X_n+1 = X_n + h*(I - hJ)^-1*F(X_n)


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

 

[ametista82]

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:

f
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
[\code]


How can I implement this with implicit integration?

thank you..

 

[matematikawan]

I'm not that familier with C++ code. If you could express your problem mathematically, probably we could come out with the algorithm.

 

[ametista82]

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?

 

[matematikawan]

I'm not sure whether I properly understand your problem here. I think your r is also a dependent variable in t. Epsilon and m are constants.

V_n+1 = V_n + dt*f(r_n+1,t_n+1)

We need another DE in dr/dt.