Solved

Posted on 1998-02-19

could anyone please help me with source code or some advice on writing a program to solve a set of 3 equations with 3 unknowns. i have found code for solving sets of 2 equations with 2 unknowns but i can't seem to convert this to deal with 3 equations.

32 Comments

A LA equation is of the form

Ax = y where A is the 3x3 coefficient matrix

x are the set of unknowns

y are constants

OK?

If A is nonsingular (positive semi definite) there is a unique solution. Other wise the solution is parametric.

I assume that you have a uniquely solvable problem.

To solve this we should find the inverse of matrix A, then the solution will be:

x = ((A)^-1).y here (A)^-1 denotes the inverse of A

and

(A)^-1 = Transpose(adj(A))/Det(A)

adj is the adjoint which I can hardly explain in plain text

tranpose is normal martix transposition

Det is normal martix determinant

SO,

typedef

float TVector[3];

TVector *PVector;

TVector TMatrix[3];

TMatrix *PMatrix;

// The A is the coefficient matrix

// y is the constant matrix that is formed by each

// equations constant

// The solution to problem returns in y where

// y = [x1 x2 x3]

int Solve3_3(PMatrix A, PVector y) {

float x1, x2, x3;

float a,b,c,d,e,f,g,h,i;

float det; //determinant

TMatrix I; //inverse matrix would be here

// now invert the matrix

// for easy programming transfer the matrix

a = A->[0][0]; b = A->[0][1]; c = A->[0][2];

d = A->[1][0]; e = A->[1][1]; f = A->[1][2];

g = A->[2][0]; h = A->[2][1]; i = A->[2][2];

// Find the determinant

det = a*e*j + d*h*c + g*b*f - (g*e*c + h*f*a + j*d*b);

// the solution does not exists if det == 0

if (det == 0) return 0;

// Find adjoint transpose this is a bit lengthy

I[0][0] = e*j - f*h; I[0][1] = b*j - c*h; I[0][2] = b*f - e*c;

I[1][0] = d*j - f*g; I[1][1] = a*j - c*g; I[1][2] = a*f - d*c;

I[2][0] = d*h - e*g; I[2][1] = a*h - b*g; I[2][2] = a*e - d*b;

// find the solution

x1 = I[0][0]*y->[0]+ I[0][1]*y->[1]+ I[0][2]*y->[2];

x2 = I[1][0]*y->[0]+ I[1][1]*y->[1]+ I[1][2]*y->[2];

x3 = I[2][0]*y->[0]+ I[2][1]*y->[1]+ I[2][2]*y->[2];

// transfer solution to return vector

// note that we do not yet divice by det so do that

y->[0] = x1 / det;

y->[1] = x2 / det;

y->[2] = x3 / det;

// return success

return 1;

};

Test the code please, I do not have C compiler on this machine. But it should be %99 correct. To convice you I am doing PhD in electronics. Bugs are wellcome, I do solve them

Please comment immediatelly

Igor

Does that help you any?

i) -1/y[ax - x^2 - bxz + (x^2)z - ab + bx + (b^2)z - axz] = 0

ii) -1/y[az - xz - bz^2 + xz^2 - a + x + bz -xz] = 0

iii) -85/2y + 1/2y^2[a^2 - 2ax + x^2 - 2abz + 2(a+b)xz -2bxz^2

+ b^2z + x^2z^2 = 0

Try looking at these softwares on your campus computer if you are a college student.

can you start by solving

[ax - x^2 - bxz + (x^2)z - ab + bx + (b^2)z - axz] = 0

[az - xz - bz^2 + xz^2 - a + x + bz -xz] = 0

method (simple matrix mathematics)..

It should be easy to implement as code.. :)

Full Fortran code for solving N nonlinear equations in N variables. Should transfer to c/c++ fairly easily, if you

watch your implicit conversions.

It actually solves the equation for

Ax=b where A,x and b are all matrix

eg.

a11x1 + a12x2+ a13x3 =b1

a21x1 + a22x2 +a23x3 =b2

a31x1 + a31x2 +a33x3 =b3

the above is a 3 by 3 matrix that you desire

where

A is a 3 by 3 matrix

x a 3 by 1 matrix

and b a 3 by 1 matrix

Eg

8x2 + 2x3 =-7

3x1 + 5x2 + 2x3 = 8

6x1 + 2x2 + 8x3 =26

We have to arrange the matrix first..always

That is to make sure that the 1st row will have a a11

the 2nd row a22. this is called pivoting

We set the pivot to x1 first

Pivoting

3x1 + 5x2 + 2x3 = 8 Pivot

8x2 + 2x3 =-7

6x1 + 2x2 + 8x3 =26

The second step is called elimination

That is start by eliminating the x1

To do this we multiply the first(pivot) row by 2 and use this

answer to subtract the 3rd row and then place this new

answer in row 3.

3x1 + 5x2 + 2x3 = 8

8x2 + 2x3 =-7

-8x2+ 4x3 =10

The second step is repeated and now we eliminate

the next variable, that is x2. The pivot is now

set to the 2nd row

This is done by adding row2 to row3 and placing the

answer in row3

3x1 + 5x2 + 2x3 = 8

8x2 + 2x3 =-7

6x3 = 3

Lastly we do back substitution to find

the rest of the answer

You would notice that the pivoting and

elimination actually only serve one

purpose, that is to obtain a triangular

matrix.

Regards

Ben

xtac@pacific.net.sg

thanks anyway.

This problem can be solved by iterative techiques using optimization.

if we have

f(x,y,z) = 0;

g(x,y,z) = 0;

h(x,y,z) = 0;

and from the equations you specified we know that y <> 0

The question can be translated into form

Minimize(f(x,y,z) + g(x,y,z) + h(x,y,z))

s.t.

z <> 0

Does this mean something,(or do I continue?)

Igor

Minimize (

Yes, you are right, make them convex and continious. But we have a second order pole at y=0. Any way we may get rid of it.

Come on lets discuss and finish this today

Igor

What descent algorithm do we use -in your opinion- ozo? Steepest descent, conjugate gradient...?

Igor

For smooth functions like these, I'd think most any descent algorithm should suffice.

There can't be that many local minima, so just start over if it gets stuck.

But is rskell looking for any (approximate) solution,

or for an equation giving all solutions?

As you said in the early comments we can safely solve only the first two equations with two unknowns then find z that (x,z) satisfying the third, if only rskell give us those a and b. Alternatively, since we eliminate y, we can find solutions for (x,z) on 2tupples in equation i; get the only ones that fits into eqn2, and jump to the eqn 3 with them. As known, in general we can solve a single eqn with 2 unknowns if it is nonlinear. Isn't it?(To simply attack the problem we may use this method. Theoreticaly one may raise several arguments but we can try.) I am now giving the partial derivatives, please check them out:(after muliplying them with y and grouping)

Eqn f=(I(x,z))^2

Ã° Ã°

---(f) = 2I(x,z)(1-z)(a+b-2x) , ---(f) = 2I(x,z)(x^2-(a+b)x + b^2)

Ã°x Ã°z

Eqn g=(II(x,z))^2

Ã° Ã°

---(g) = 2II(x,z)(-2x + 1) , ---(g) = 2II(x,z)(a + (b-x)(1-2z))

Ã°x Ã°z

We may start the either one, are they ok?

Igor

One is Newton-Raphson Method

The other is Globally convergent Methods.

Has anyone heard about these two Methods?

Yes, the Newton-Raphson(NR) method is a well known method for finding roots of roots of equations. It converges to the solution quite fast. It utilizes 2nd order derivatives together with the tangent hyperplane intersection with the hyperplanes formed by the axises. But when the problem is multi dimensional, we can not directly apply the 1D NR. We should form a matrix called Hessian which contitues the 2nd order partial derivatives of all functions w.r.t. all the variables in the equation set. All this are examined in the scope of Optimization(minimization or maximization of functions) and Optimum Control (minimum time, minimum energy etc...problems). The other thing "Globally convergent Methods" are include several methods. They are proposed mainly for complex error surfaces which have several local minima (for minimization problems). The well know Globally convergent method is e.g. Simulated Annealing, which gives small disturbance to the solution when a local minima is reached hoping to step over the nearest local maxima and reach a more appropriate local minima, by this way it seek the global maxima. The other are Evolutionary approaches, Boltzman machine approach and combination of them such as GESA(Guided Evolutionary Simulated Annealing).

Sorry, if I go out the scope of the question

Bye

Igor

Title | # Comments | Views | Activity |
---|---|---|---|

isEverywhere challenge | 19 | 49 | |

either24 challenge | 19 | 73 | |

static class | 3 | 44 | |

Eclipse IDE - Cannot copy/paste from console output | 8 | 15 |

This article will show, step by step, how to integrate R code into a R Sweave document

The viewer will learn how to clear a vector as well as how to detect empty vectors in C++.

Join the community of 500,000 technology professionals and ask your questions.

Connect with top rated Experts

**23** Experts available now in Live!