Numerical Methods, need a 3D Jacobian

[Feodalherren]

Homework Statement

I'm trying to write a program to solve a system of 3 non-linear equations using the Newton-Raphson method. I'm stuck on trying to figure out what the formula for a system of 3 unknowns is. I can't remember the derivation at all and after endless hours of googling and looking at pdfs (they all list a two-by-two ARG!)

Homework Equations

The Attempt at a Solution

Does anyone know what the formula looks like for a 3by3 or an arbitrary system?
Thanks!

 

[SteamKing]

Homework Statement

I'm trying to write a program to solve a system of 3 non-linear equations using the Newton-Raphson method. I'm stuck on trying to figure out what the formula for a system of 3 unknowns is. I can't remember the derivation at all and after endless hours of googling and looking at pdfs (they all list a two-by-two ARG!)

Homework Equations

The Attempt at a Solution

Does anyone know what the formula looks like for a 3by3 or an arbitrary system?
Thanks!

The Jacobian matrix is well-defined:

https://en.wikipedia.org/wiki/Jacobian_matrix_and_determinant

I assume you have a system of 3 non-linear equations in 3 unknowns.

Call the unknowns x₁, x₂, and x₃. There will be three functions f(x₁, x₂, x₃). Call the three functions f₁, f₂, and f₃.

This article shows how to set up the N-R method (see p. 9; set n = 3):

http://ocw.usu.edu/Civil_and_Enviro...ivil_Engineering/NonLinearEquationsMatlab.pdf

You will have to find the inverse J^-1 of the Jacobian, or better still, put your system of non-linear N-R equations in matrix form Ax = B and solve this system to avoid having to compute J^-1

 

[Chestermiller]

You will have to find the inverse J^-1 of the Jacobian, or better still, put your system of non-linear N-R equations in matrix form Ax = B and solve this system to avoid having to compute J^-1

Isn't solving that equation the same as inverting the Jacobian, since A is the Jacobian?

Chet

 

[SteamKing]

Isn't solving that equation the same as inverting the Jacobian, since A is the Jacobian?

Chet

Algebraically, it is, but computing a matrix inverse numerically is more complicated and leads to large round off errors in certain circumstances. From a numerical standpoint, it is recommended that computing a matrix inverse be avoided as much as possible, and instead, the problem reformulated to solve a system like Ax = B rather than compute A^-1B = x.

 

[Chestermiller]

Algebraically, it is, but computing a matrix inverse numerically is more complicated and leads to large round off errors in certain circumstances. From a numerical standpoint, it is recommended that computing a matrix inverse be avoided as much as possible, and instead, the problem reformulated to solve a system like Ax = B rather than compute A^-1B = x.

When you refer to getting the matrix inverse, I think you're talking about using determinants to get the inverse. Of course, there are other ways of getting the inverse, such as solving Ax=B using ordinary elimination techniques. If the person really had to know the inverse of A, they could simultaneously apply the same sequence of elimination operations to the identity matrix I.

Chet

 

[SteamKing]

When you refer to getting the matrix inverse, I think you're talking about using determinants to get the inverse. Of course, there are other ways of getting the inverse, such as solving Ax=B using ordinary elimination techniques. If the person really had to know the inverse of A, they could simultaneously apply the same sequence of elimination operations to the identity matrix I.

Chet

Even calculating the inverse of a matrix using elimination is advised against unless absolutely necessary. Roundoff may not become significant in all cases, but the number of operations required to calculate the inverse grows proportional to n³, where n is the order of the matrix. In solving large problems involving hundreds, if not thousands of equations, no one forms the inverse; the solutions are typically obtained using some form of elimination or iteration.

Calculating the inverse of a matrix involves solving the system AA^-1 = I, where I is the n × n identity matrix, A is the n × n matrix to be inverted, and A^-1 is the unknown n × n inverse of A.

Most of the books on numerical analysis I've read, including Forsythe, Malcolm, and Moler, Computer Methods for Mathematical Computations, advise against computing the inverse unless it is unavoidable.

http://www.netlib.org/lapack/lawnspdf/lawn27.pdf {See p. 2}

In the N-R method, the Jacobian matrix is updated as the solution vector is updated during the iterations. Since the inverse of the Jacobian is required only as a means of calculating a new solution vector, why not just solve the equations to obtain the new solution vector, instead of slavishly following the matrix equation?

It's why we use the quadratic formula to solve a quadratic equation, instead of laboriously completing the square for every equation. In the end, we get the same result, but one method involves fewer steps in the calculation.

 

[Feodalherren]

Thanks a lot guys! One follow up question: what exactly do you mean when you're saying set it up as Ax=B. I'm not sure how I'd do that with a function that involves transcendental equations.

Can you give me a concrete example with some random non-linear function of, say 2 variables?

Thanks again!

Edit: I still can't figure out how they got the numerator. That pdf explains how to get the denominator clearly but the numerator still seems completely arbitrary.

 

[Chestermiller]

If f is the vector of non-linear functions you are trying to get to zero and x_n is the approximation to the solution vector after n applications of NR, then expanding f in a Taylor series about x_n, we obtain:

f(x) = f (x_n)+J(x_n)(x - x_n)

We are trying to make f(x_n+1) equal to zero. So we write:

f (x_n)+J(x_n)(x_n+1 - x_n) = 0

or
J(x_n)(x_n+1 - x_n)=-f (x_n)

So, A = J(x_n), B= -f (x_n) and the unknown is (x_n+1 - x_n).

Hope this makes sense.

Chet

 

[Feodalherren]

Now I'm even more confused, how can you have x_n+1 in the equation if that's what I'm trying to find? My code so far looks like this:

MATLAB:

while count<100 && ea>es
    count=count+1;

% The Jacobian Matrix
J=[df1x(x) df1y df1z ; df2x df2y(y) df2z ; df3x df3y df3z(z)];
D=det(J); %the denominator in Newton-Raphson is the determinant of the Jacobian

%use the Newton-Raphson formula to find new roots
newx = x - ((f1(x,y,z)*df2y(y) - f2(x,y,z)*df1y)/D);
newy = y - ((f2(x,y,z)*df1x(x) - f1(x,y,z)*df2x)/D);
newz = z - ((f1(x,y,z)*df2y(y) - f2(x,y,z)*df1y)/D);
%compute error for xyz
eax=abs((newx-x)/newx)*100;
eay=abs((newy-y)/newy)*100;
eaz=abs((newz-z)/newz)*100;
A=[eax eay eaz];

ea=max(A); %error

x=newx;y=newy;z=newz; %re-define for next iteration
end %while

I know that the newx,newy and newz terms are wrong.

 

[Chestermiller]

Now I'm even more confused, how can you have x_n+1 in the equation if that's what I'm trying to find? My code so far looks like this:

MATLAB:

while count<100 && ea>es
    count=count+1;

% The Jacobian Matrix
J=[df1x(x) df1y df1z ; df2x df2y(y) df2z ; df3x df3y df3z(z)];
D=det(J); %the denominator in Newton-Raphson is the determinant of the Jacobian

%use the Newton-Raphson formula to find new roots
newx = x - ((f1(x,y,z)*df2y(y) - f2(x,y,z)*df1y)/D);
newy = y - ((f2(x,y,z)*df1x(x) - f1(x,y,z)*df2x)/D);
newz = z - ((f1(x,y,z)*df2y(y) - f2(x,y,z)*df1y)/D);
%compute error for xyz
eax=abs((newx-x)/newx)*100;
eay=abs((newy-y)/newy)*100;
eaz=abs((newz-z)/newz)*100;
A=[eax eay eaz];

ea=max(A); %error

x=newx;y=newy;z=newz; %re-define for next iteration
end %while

I know that the newx,newy and newz terms are wrong.

You are going to be solving the last (linear algebraic) equation I wrote for x_n+1.

If you are insistent on using determinants to invert J, then you can set the calculation up manually to do that. Or you can have a software subroutine to invert J. Or you can use a linear equation solver to solve the corrector equation for x_n+1. SteamKing and I both recommend that you don't use determinants to invert J. It is much preferred that you use a linear equation solver, which employs techniques that result in much less roundoff error.

Chet

 

[Feodalherren]

So then
x_n+1=x_n - f(x_n / J(x_n))
?

which I could code as:

newx=x-f(x,y,z)/J(x,y,z)

But then what was the whole discussion about not calculating J^-1?

Thanks for the help! It's much appreciated.

 

[SteamKing]

So then
x_n+1=x_n - f(x_n / J(x_n))
?

which I could code as:

newx=x-f(x,y,z)/J(x,y,z)

But then what was the whole discussion about not calculating J^-1?

Thanks for the help! It's much appreciated.

Your equation only works for one non-linear equation in several unknowns. It's essentially a re-statement of Newton's method.

When you have a system of n non-linear equations in n unknowns, finding the new values of (x₁, x₂, ..., x_n) requires a different technique, one involving finding the simultaneous solution of n equations for the n unknowns.

 

[Chestermiller]

Your equation only works for one non-linear equation in several unknowns. It's essentially a re-statement of Newton's method.

When you have a system of n non-linear equations in n unknowns, finding the new values of (x₁, x₂, ..., x_n) requires a different technique, one involving finding the simultaneous solution of n equations for the n unknowns.

Actually, his equation looks OK to me if, by 1/J he means J^-1.

Chet

 

[SteamKing]

Actually, his equation looks OK to me if, by 1/J he means J^-1.

Chet

Maybe, but I got the feeling that the OP may have missed the distinction that 1/ J and J^-1 represent two different things.

 

[Feodalherren]

Your equation only works for one non-linear equation in several unknowns. It's essentially a re-statement of Newton's method.

When you have a system of n non-linear equations in n unknowns, finding the new values of (x₁, x₂, ..., x_n) requires a different technique, one involving finding the simultaneous solution of n equations for the n unknowns.

That's exactly what I'm trying to do, but I can't find the formula for it. This is beginning to vex me. This class isn't a math class why does the professor not give us formulas?! I don't have time to go back and re-learn linear algebra for a programming class just so that I can derive a stupid equation.
Does nobody have a link or something to a website that will just give me the equation that I need?

Oh and yeah I did mean the inverse of the Jacobian, as I said I'm rusty on my linear algebra and that includes the notation.

 

[SteamKing]

That's exactly what I'm trying to do, but I can't find the formula for it. This is beginning to vex me. This class isn't a math class why does the professor not give us formulas?! I don't have time to go back and re-learn linear algebra for a programming class just so that I can derive a stupid equation.
Does nobody have a link or something to a website that will just give me the equation that I need?

Oh and yeah I did mean the inverse of the Jacobian, as I said I'm rusty on my linear algebra and that includes the notation.

Welcome to the Real World. Often, you have to do a little research to find the answers you need. Sometimes, you even have reach back to classes finished long ago for the necessary information and techniques.

---------------------------------------------

It's not clear what you mean here:

That's exactly what I'm trying to do, but I can't find the formula for it.

Are you talking about how to solve an n × n system of equations? Since you are dealing with a 3 × 3 system, Cramer's Rule is often used for hand calculations. However, given the iterative nature of the N-R method, cranking thru all those 3 × 3 determinants by hand may be even more vexing, given your current state of mind.

https://en.wikipedia.org/wiki/Cramer's_rule

The Rule of Sarrus gives the nuts and bolts of evaluating 3 × 3 determinants.

https://en.wikipedia.org/wiki/Rule_of_Sarrus

This article:

http://ocw.usu.edu/Civil_and_Enviro...ivil_Engineering/NonLinearEquationsMatlab.pdf

also contains some Matlab routines which could be used to find the solutions to your system of equations. IDK if you have access to Matlab, but it's a place to start. See pp. 10 and following for the routines. You may have to adapt what is shown for your particular 3 × 3 system of equations.

 

[Chestermiller]

SteamKing and I have been pleading with you to consider using another method than determinants to solve your set of 3 linear algebraic equations in 3 unknowns. Don't you have any linear equation solves available to you in your software library? All you need to do is to feed the Jacobian matrix and the right hand side to the linear equation solver, and it spits out the solution vector. It's a zillion times better than writing your own code to do the same thing.

Chet