#### Archived

This topic is now archived and is closed to further replies.

# Implicit Euler

This topic is 5554 days old which is more than the 365 day threshold we allow for new replies. Please post a new topic.

## Recommended Posts

Hi, I''m trying to learn about Implicit Euler integration so I can implement it in my own physics engine. I have found some documents that mention it but they never go as far as giving details of what the equations are or how to implement it. I was wondering if someone could point me to some docs or books or even give a brief explanation of what it is and how I could go about implementing it. I am familiar with Euler (explicit), mid-point, Runge-Kutta and verlet integration methods. Thanks for your help, lusinho

##### Share on other sites
I'm not sure what you mean by "implicit Euler" integration. Implict formulae are those like xy = 1, x^2 + y^2 = 2, etc. In general they're of limited interest to games programmers as they are difficult and expensive to deal with in code.

As for Euler, he was one of the most significant contributors to modern maths, particularly calculus, and his name is associated with many interesting discoveries. A quick seach of MAthworld turned up the Euler-Lagrange Differential Equation, the Euler Integral, the Euler Differential Equation and many more.

[edited by - johnb on November 21, 2002 10:30:34 AM]

[edited by - johnb on November 21, 2002 10:30:46 AM]

##### Share on other sites
If you can find the David Baraff tutorials (which I''m sure google can) I think he had some implicit euler stuff. And I think there where in the later sections...but I could be wrong on both scores.

Good luck

##### Share on other sites
http://www.physics.udel.edu/faculty/macdonald/Ordinary%20Differential%20Equations/The%20Implicit%20Euler%20Method.htm

##### Share on other sites

Yes, you can find the latest version of Baraff''s tutorials at

http://www.pixar.com/companyinfo/research/pbm2001/index.html

and if you look under Implicit Methods you''ll find that these come about to solve stiff ODEs.

Another name for Implicit Euler is Backward Euler (I think).

As I understand it the normal (explicit) Euler integration method can be expressed like

X1 = X0 + h * Y(X0)

where X1 is the new position (orientation), X0 the current position, h is the step (delta time) and Y(X0) is the first derivative of X0 (velocity).

For Implicit Euler the change is simple

X1 = X0 + h * Y(X1)

so the only difference is that the derivative is taken on the new state and not the current one.
But because we don''t know the new state yet, then you can''t use this function directly and have to do it another way instead.

I can''t really follow the whole explanation from Baraff, and I have a feeling that even if I did that still wouldn''t help me to implement this algorithm.
This is becasue in some other documents I''ve read that you need to calculate some big matrices and other complicated things in order to get to the final result.
I can''t actually remember which documents these were, but I do remember them only mentioning it but never explaining how to do it.

##### Share on other sites
>so the only difference is that the derivative
>is taken on the new state and not the current one.
>But because we don't know the new state yet, then you
>can't use this function directly and have to do it another way >instead.

You're correct. Let's write the equation this way:

(I denote X0 with X)

X1 = X + h*Y(X1) // replace X1 with X+DX
X+DX = X + h*Y(X+DX) // subtract X
DX = h*Y(X+DX)

the problem is that we are not able to compute Y(X+DX) directly. Therefore, we use a technique called linearization, that is, we approximate Y(X+DX) using the normal Euler method:

Y(X+DX) ~= Y(X) + (@Y/@X)*DX // ~= is "approximates"

Here @Y/@X is the partial derivative of Y wrt. X. If Y is a vector ("forces of particles") and X is a vector ("positions of particles"), the partial derivative will be a matrix. The most complicated part of the algorithm is computing this matrix.

Finally we can solve for DX:

DX = h*Y(X) + h*(@Y/@X)*DX // solving for DX

[I-h*(@Y/@X)]*DX = h*Y(X) // I denotes the identity matrix

If X is a vector, and Y is a vector, this becomes a linear system of equations. Usually people use the so called conjugate gradient method for solving the equation. After it's done you have found DX: the problem of backward Euler has been successfully solved!

>This is becasue in some other documents I've read that you need >to calculate some big matrices and other complicated things in >order to get to the final result.

True: the matrix (@Y/@X) is a "big" matrix. If the size of X and Y (both assumed to vectors) is n, the matrix will be of size nxn. The lucky thing is, that in the matrix, most nodes are zero, so you can save a hell lot of memory and computational power. The matrix @Y/@X is defined as follows: let

Y = (y1, y2, y3, ..., yn)
X = (x1, x2, x3, ..., xn)

Then the matrix is defined as follows:

@Y/@X = (@Y/@x1, @Y/@x2, @Y/@x3, ..., @Y/@xn)

(where each node @Y/@xi is a column vector so @Y/@X really is a matrix)

More generally, we can write out the (j,i) node of the matrix as:

(@Y/@X)ji = @fj/@xi

Hope this rant helps a bit.

- Mikko Kauppila

[edited by - uutee on November 21, 2002 4:27:53 PM]

##### Share on other sites

Thanks Mikko Kauppila, I think that helps as now I can start learning about the Conjugate Gradient Method and how to implement that.

Doing a quick seach on the web I found that one of Baraff''s papers actually talks about this too. It''s "Large steps in cloth simulation" at http://www.pixar.com/companyinfo/research/deb/sig98.pdf

One of his refenrences is "An introduction to the Conjugate Gradient Method without the agonizing pain" at

I''ll need some time to digest both of these documents.

Thanks,
lusinho

##### Share on other sites
Baraff''s paper (Large steps) explains the subject nicely for cloth simulation. There''s also a more general paper on the subject by Baraff, it''s in Siggraph course notes. Google for "physically based modeling baraff", you should find an article called "implicit methods for differential equations" quickly. "Large steps in cloth simulation" explains implicit integration in a great, great fashion, but it explains the cloth forces quite badly, to speak the truth.

Shewchuk''s paper on the CG method is truely a great piece of information. It gets quite hard to read after a couple of tens of pages though (as it approaches the CG method), but the "Canned algorithms" section of the paper gives a nice pseudo code of the CG algorithm, if you just need to implement it. Roughly speaking CG solves the equation A*x=b (where A is a matrix, b is a vector and x is the vector we''re trying to solve) by starting with an initial guess vector x0 (for example, the null vector) and by iteratively adjusting it. Each iteration requires one multiplication of a vector by A, which is a quick-to-do operation for sparse matrices (such as the one''s we tend to have in implicit methods).

- Mikko Kauppila

##### Share on other sites
quote:
Original post by lusinho
Another name for Implicit Euler is Backward Euler (I think).

Yep, you''re right about that.

Graham Rhodes
Senior Scientist
Applied Research Associates, Inc.

##### Share on other sites
After reading those documents I am slowly crawling towards a better understanding and I have a few more questions and coments that I would appreciate any views.

From Baraff''s papers our state vector is

	X(t) =	[ x(t) ]	position		[ q(t) ]	orientation		[ P(t) ]	linear momentum		[ L(t) ]	angular momentum

The derivative of the state vector is

	d         .		[ v(t) ]	linear velocity	-- X(t) = X(t) =	  .	dt			[ q(t) ]	angular velocity				[ F(t) ]	force				[ T(t) ]	torque

The final Implicit Euler equation is

X(t+h) = X(t) + delta_X

where delta_X is given by

	  1       d   .                    .	[ - * I - --( X(t) ) ] * delta_X = X(t)	  h       dX

and as this equation is of the form

A * x = b

where b is a known vector, A is a known matrix (square, symmetric and positive-definite) and x is our unknown vector, we can solve it using the Conjugate Gradient Method.

Great, but I don''t know what the A matrix is.

Baraff tells us that A is a nxn sparse matrix where n is the size of the state vector and also that the (i,j) entry only exists iff (if and only if)

f(i) depends on X(j)

After some head-scratching I thought I understood what this meant and made an attempt at building this matrix.

A =[   1 	d( v(t) )	d( v(t) )	d( v(t) )	d( v(t) )	][   - -	---------	---------	---------	---------	][   h	d( x(t) )	d( q(t) )	d( P(t) )	d( L(t) )	][									][	   .		   .		   .		   .		][ 	d( q(t) )   1	d( q(t) )	d( q(t) )	d( q(t) )	][	---------   - -	---------	---------	---------	][	d( x(t) )   h	d( q(t) )	d( P(t) )	d( L(t) )	][									][ 	d( F(t) )	d( F(t) )   1	d( F(t) )	d( F(t) )	][	---------	---------   - -	---------	---------	][	d( x(t) )	d( q(t) )   h	d( P(t) )	d( L(t) )	][									][ 	d( T(t) )	d( T(t) )	d( T(t) )   1	d( T(t) )	][	---------	---------	---------   - -	---------	][	d( x(t) )	d( q(t) )	d( P(t) )   h	d( L(t) )	]

Then I tried to identify which elements would be zero according to the above rule.
For each element I thought "does the term on top depend on the term on the bottom?".
So for example for the element (0,0): "does velocity depend on position?"
The answer is no. Position depends on velocity but not the other way around.

So this is the matrix I arrived to

A =[	1 		d( v(t) )			][	-	0	---------	    0		][	h		d( P(t) )			][							][	   	   	   		   .		][ 		1			d( q(t) )	][	0	-	    0		---------	][		h			d( L(t) )	][							][ 			    1				][	0	0	    -		    0		][			    h				][							][ 					    1		][	0	0	    0		    -		][					    h		]

This looks a bit too sparse.

Is my initial matrix correct?
Did I make a mistake in identifying the zero components?
How do I calculate the remaining elements of the matrix:

d( v(t) )--------- = ?d( P(t) )   .d( q(t) )--------- = ?d( L(t) )

Thanks for any help.
lusinho

##### Share on other sites
It seems you''re trying to implement implicit integration for rigid body simulation. The thing is, that this is not very useful, because the forces rigid bodies get are rarely dependent on the other rigid bodies. Therefore, as there''s very little interaction between the bodies, the matrix ends up really sparse.

The backward Euler method is useful on systems where we have a great number of objects (usually particles) which are connected somehow (for example springs). Cloth simulation is a great example of this. The matrix will still end up being sparse (on the average 8 nodes on each row).

You derived your matrix correctly. You''re correct, that kind of "does the quantity we''re derivating depend on the quantity with which we''re derivating with respect to?" is an intuitive way to understand the sparsity of the generated matrices.

I don''t know what you''re thinking, but the implicit Euler method is definitely NOT to be applied on the bodies one by one. The whole idea of the method is that it supports GLOBAL interactions between the objects. Therefore we need to construct one, huge but luckily sparse matrix.

- Mikko Kauppila

##### Share on other sites
Well, my idea was to understand the whole process for a single rigid object and then extend it for more, but I wasn''t aware that in this case you need to take into account all interacting objects.

I planned to use this method for single rigid bodies, multiple rigid bodies linked by joints (articulated), cloth and maybe other things I find I can do with this method.

The reason why I decided to explore this method is because my current integration methods fail under certain circumstances, especially under the influence of springs.

I looked briefly into LCP based methods and read another of Baraff''s papers ("Linear Time Dynamics using Lagrange Multipliers") and it just went flying straight above my head so I thought I would tackle something else.

Ok, I''m not going to give up on this one though.

Let''s take the example of cloth and consider three particles then (I assume this could be a minimum number so I can understand the process without making it too complicated).

So because we are considering particles then our state vector is reduced as there is no more orientation component:

	X(t) =	[ x(t) ]	position		[ P(t) ]	linear momentum

The derivative of the state vector is

	.	X(t) =	[ v(t) ]	linear velocity		[ F(t) ]	force

If we consider three particles such as particle A is linked to B, and B to C, but A is not linked to C

	A----B----C

We can put all the state vectors into a joint state vector and the same for the derivative vector

				.	X(t) =	[ A.x(t) ]	X(t) =	[ A.v(t) ]		[ A.P(t) ]		[ A.F(t) ]		[ B.x(t) ]		[ B.v(t) ]		[ B.P(t) ]		[ B.F(t) ]		[ C.x(t) ]		[ C.v(t) ]		[ C.P(t) ]		[ C.F(t) ]

The Big Matrix (which I''ll call BigM - was A previously) will now be like a grid where all the particles cross.

I''m just thinking out loud here so I may be completely wrong in which case I''d appreciate if you could correct me.

BigM = [	AA	AB	AC	][	BA	BB	BC	][	CA	CB	CC	]

where the particles are not interacting then those components will be zero.

BigM = [	AA	AB	0	][	BA	BB	BC	][	0	CB	CC	]

For the diagonal components then we have the previous A matrix from the previous post.
For example the AA matrix would be:

AA =[	1 	d( A.v(t) )	][	-	-----------	][	h	d( A.P(t) )	][				][ 		    1		][	0	    -		][		    h		]

For the components where the particles interact then we should have something like this.

AB = [   1 	d( A.v(t) )	d( A.v(t) )	][   - -	-----------	-----------	][   h	d( B.x(t) )	d( B.P(t) )	][					][	d( A.F(t) )	1   d( A.F(t) )	][	-----------	- - -----------	][	d( B.x(t) )	h   d( B.P(t) )	]

So at this point we could think "is the velocity of A dependant on the position of B?" and maybe some of these would be quite obvious but maybe that depends on the type of relationship that the particles have.
I remember reading that there is different types of springs for the cloth models and maybe depending on the type then this matrix could be different (just guessing here).

Whether I''m right or not I will still need to calculate expressions like:

d( A.v(t) )----------- = ?d( A.P(t) )d( A.v(t) )----------- = ?d( B.P(t) )

I have a feeling that these may have to be implemented as a case by case, because as I said above they may be different for the different types of spring (I think).
There must be documentation that explains this, or maybe Baraff''s paper explains it and I just didn''t get it. I''ll have to read it again.

Thanks again Mikko.
lusinho

##### Share on other sites
>There must be documentation that explains this,
>or maybe Baraff''s paper explains it and I just
>didn''t get it. I''ll have to read it again.

Baraff''s paper doesn''t explain construction of the matrix very well (quite badly to speak the truth).

>X(t) =
>[ x(t) ]position
>[ P(t) ]linear momentum

It''s of course possible to derive the equations this way by constructing a state vector which is a combination of positions and linear momentums. I haven''t implemented it myself. It''s slower and more complex. We get off MUCH more easily if we state in the original implicit Euler formulation:

deltaX = (v+deltaV)*timeStep
deltaV = f(x+deltaX,y+deltaY)

From this we only solve deltaV. deltaX is then trivial to compute. Computing deltaV only results in a simpler and smaller matrix. This is explained really well in Baraff''s paper (Large Steps [98] and the one from "physically based modeling principles and practice".

>where the particles are not interacting then those components >will be zero.
>
>BigM =
>[AA AB 00]
>[BA BB BC]
>[00 CB CC]

This is not just a bit but completely right! You''ve understood the method this far perfecly. In cloth simulation the final matrix is a sum of three matrices, A=M-h*dFdV-h^2*dFdX (they all have the same sparsity patterns, so summing them is easy). Let''s concentrate on the dFdX matrix (force derived wrt. position) which is the only "compulsory" matrix (dFdV is not compulsory if you have no damping forces in the first place, which I strongly recommend). In this case, AA in BigM means "force exerted on A derived wrt. position of A" and so on.

>d( A.v(t) )
>----------- = ?
>d( A.P(t) )

The actual derivatives are the hardest ones to compute. By solving only deltaV all we need to compute is derivatives of forces wrt. positions (and with damping forces forces wrt. velocities). External forces (mouse interactions) don''t affect the derivatives. Only if the force depends on the positions/velocities of the particles do they affect the derivatives. The only forces of this kind are the spring forces (I personally handle air friction by multiplying velocities with 0.999 at the end of the step

So in the end, the only problem left is to derivate the spring forces. The formulas are simple, and it''s shown in this paper: Kwang-Jin Choi, Hyeong-Seok Ko: Stable but responsive cloth. ACM Transactions on Graphics 21(3): 604-611 (2002)

- Mikko Kauppila