1. Limited time only! Sign up for a free 30min personal tutor trial with Chegg Tutors
    Dismiss Notice
Dismiss Notice
Join Physics Forums Today!
The friendliest, high quality science and math community on the planet! Everyone who loves science is here!

Numerical Methods, need a 3D Jacobian

  1. Oct 7, 2015 #1
    1. The problem statement, all variables and given/known data
    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!)

    2. Relevant equations
    Untitled.jpg

    3. The attempt at a solution

    Does anyone know what the formula looks like for a 3by3 or an arbitrary system?
    Thanks!
     
  2. jcsd
  3. Oct 7, 2015 #2

    SteamKing

    User Avatar
    Staff Emeritus
    Science Advisor
    Homework Helper

    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 x1, x2, and x3. There will be three functions f(x1, x2, x3). Call the three functions f1, f2, and f3.

    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
     
  4. Oct 7, 2015 #3
    Isn't solving that equation the same as inverting the Jacobian, since A is the Jacobian?

    Chet
     
  5. Oct 7, 2015 #4

    SteamKing

    User Avatar
    Staff Emeritus
    Science Advisor
    Homework Helper

    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.
     
  6. Oct 8, 2015 #5
    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
     
  7. Oct 8, 2015 #6

    SteamKing

    User Avatar
    Staff Emeritus
    Science Advisor
    Homework Helper

    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 n3, 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.
     
  8. Oct 8, 2015 #7
    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.
     
    Last edited: Oct 8, 2015
  9. Oct 8, 2015 #8
    If f is the vector of non-linear functions you are trying to get to zero and xn is the approximation to the solution vector after n applications of NR, then expanding f in a Taylor series about xn, we obtain:

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

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

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

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

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

    Hope this makes sense.

    Chet
     
  10. Oct 8, 2015 #9
    Now I'm even more confused, how can you have xn+1 in the equation if that's what I'm trying to find? My code so far looks like this:

    MATLAB:
    Code (Text):

    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.
     
  11. Oct 8, 2015 #10
    You are going to be solving the last (linear algebraic) equation I wrote for xn+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 xn+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
     
    Last edited: Oct 8, 2015
  12. Oct 8, 2015 #11
    So then
    xn+1=xn - f(xn / J(xn))
    ?

    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.
     
  13. Oct 8, 2015 #12

    SteamKing

    User Avatar
    Staff Emeritus
    Science Advisor
    Homework Helper

    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 (x1, x2, ..., xn) requires a different technique, one involving finding the simultaneous solution of n equations for the n unknowns.
     
  14. Oct 8, 2015 #13
    Actually, his equation looks OK to me if, by 1/J he means J-1.

    Chet
     
  15. Oct 8, 2015 #14

    SteamKing

    User Avatar
    Staff Emeritus
    Science Advisor
    Homework Helper

    Maybe, but I got the feeling that the OP may have missed the distinction that 1/ J and J-1 represent two different things.
     
  16. Oct 8, 2015 #15
    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.
     
  17. Oct 8, 2015 #16

    SteamKing

    User Avatar
    Staff Emeritus
    Science Advisor
    Homework Helper

    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:
    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.
     
  18. Oct 8, 2015 #17
    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
     
Know someone interested in this topic? Share this thread via Reddit, Google+, Twitter, or Facebook

Have something to add?
Draft saved Draft deleted



Similar Discussions: Numerical Methods, need a 3D Jacobian
Loading...