Given two systems of linear equations with the same coefficient matrix A: A * x1 = b1; // Eq: (1) A * x2 = b2; // Eq: (2) Adding both equations yields: A * ( x1 + x2 ) = b1 + b2 // Eq: (3) Now I substitute x = x1 + x2 b = b1 + b2 and solve A * x = b // Eq: (4) Q: If I now have the solution of x, is it possible to find the solutions of individual x1 and x2 *without* solving either equation (1) or (2)? Or in other words: Does a coefficient vector alpha exist such that: x1 = alpha * x x2 = (1 - alpha) * x Would it be possible to derive the relative scale of each element of x1 and x2 from the relation of b1 to b2, e.g.: a1 = b1 / (b1 + b2) a2 = b2 / (b1 + b2) Regards, Dirk

If you want to efficiently solve multiple systems with the same coefficient matrix, you should solve for the inverse of the matrix. The relations you are describing must exist, but I think it would take more work to solve for them than it would to invert the matrix and solve each system separately.

If you know what Gauss-Jordan elimination is you can use that to solve multiple systems at the same time with the same coefficient matrix with less overhead than doing each equation separately.

When doing Gauss-Jordan elimination you perform the row operations on the coefficient matrix to get it into reduced row echelon form and while doing this you perform the same row operations on the 'b' column vector. If you have more than one system you just perform the operations on all the column vectors at the same time. It is exactly like how you would find the inverse of the matrix. However instead of augmenting the identity matrix with the coefficient matrix you augment each of the 'b' column vectors.

I had the following idea. When solving the DAE of a rigid body system you have the following system of equations:

(1) M * dv / dt = f_c + f_ext
(2) J * v = 0
(3) f_c = JT * lambda

Using a velocity constrained formulation you typically solve:

(4) J * W * JT * lambda + J * ( v(t) / dt + W * f_ext ) = 0

In the presence of contacts or other unilateral constraints you get a mixed LCP. In this form you don't counter any numerical drift. The typical solution is to add some kind of penalty (Baumgarte) term - something like c = 0.2 * error / dt and plug this term into equation (4). So finally one arrive at:

(5) J * W * JT * lambda + J * ( v(t) / dt + W * f_ext ) - c / dt = 0

This is the formulation used in the ODE. The problem of this form is that we add energy to the persistent momentum of the bodies, because of the presence of the penalty term. While this is not a problem for small timesteps and small interpenetrations, for me this is quite some problem, since I don't have nice small timesteps nor small interpenetrations and I get quite some unpretty jittering...

My first idea now was that solving equation (5) is the same like solving the following system - simple super-positioning:

(4) J * W * JT * lambda1 + J * ( v(t) / dt + W * f_ext ) = 0
(6) J * W * JT * lambda2 - c / dt = 0
(7) lambda = lambda1 + lambda2

Since I don't want to solve two MLCP I wondered if I could find lambda1 and lambda2 from the relation of the "b" terms of the system. Since I solve the LCP numerically (PGS) I actually don't compute the system matrix nor do actually build the inverse.

Q1:
If I can't compute lambda1 and lambda2 from solving only equation 5 - is it possible to solve the PGS such that I could find lambda1 and lambda2 through modifying the PGS? This way I could formulate the position update such that the error correction does NOT contribute to the persistent momentum of the bodies. This should make the constraints also a bit harder.

Q2:
There is a paper of M. Anitescu that deals with this problem as well. Does anybody know how this issue is adressed there and/or has this paper and can provide a link?


-Dirk