Hi, I am solving a sistem of lineal ecuations with the subroutine ''gaussj''. I get the answer (the inverse matrix and the solution vector) but I find not so good presicion in it.
E.g.> I can do a try with a matrix wich its inverse has somewhere 2 as exact solution, but instead the program gives me 1.999768...
I would be really happy if it tells me it is 2.000000000! :)
Does someone know what's going on and how can I improve the result?

I would be really happy if it tells me it is 2.000000000

I look at this way:

If the matrix is non-singular (that is, it has an inverse), and if the Gauss-Jordan scheme is implemented correctly, then the only way that you could get the "wrong" answer would be if roundoff error had raised its ugly head. I mean, the possibility of roundoff error is "always" there, but sometimes it is of such a low magnitude that it doesn't interfere with our expectations.

Since you apparently know what the inverse is, then the matrix is non-singular. If you are using one of the Numerical Recipes gaussj functions (and there aren't any transcription errors, and if your compiler is working correctly, etc.), then the algorithm is implemented correctly.

That leaves roundoff error.

Some matrices are more susceptible to roundoff error than others. This is discussed several times in sections of Chapter 2 in the various Numerical Recipes texts, and there is no simple answer that tells us the "best" way to deal with it for all matrices. (Well, maybe there are some simple answers, but not correct answers; not for all matrices.)

So the questions for you are:


1. Does your test program work for some systems (that is: does it give satisfactory results on the printout) but not for others?

2. What is the matrix to which you are applying gaussj? Where did it come from? Was it specially chosen as an example because it has a large condition number? Or did you just get lucky?

3. What version of Numerical Recipes gaussj are you using? Version 3? If it's not version 3, then what language? Can you share the matrix with us? How about the test program?


Regards,

Dave

[QUOTE=davekw7x;3009]I look at this way:

Hi Dave! First of all thanks so much for your replay!
Well I shall explain you a bit what it is all about...
I am building a numerical model about wings flow using a finite vortex lattice. I was told it was better to use Fortran instead of matlab cause it can take it better. So here I am.
To find a solution I have to solve a big matrix (around 70x70 to begin). But it is of course not any sort matrix. Because of the geometry of the problem it is always non-singular (all lines are lineal independent and so).

I made many test with the program I wrote and I always get a solution. I know from the problem that some numbers in the solution must be the same but instead I find little diference between them (around the 4th decimal). So I decided to do a test with a matrix which its inverse I knew beforehand. And voi la! The same problem again. So that turns all the attention in the program. I am 99% sure it is not a problem related to the matrix.

I know how gauss jordan method works, but I have no idea what the gaussj subroutine is exactly doing. As I found it in the recipie I just copied and pasted because I didnt want to program something that was made. Magicly no error was found but then I noticed this problem in the result and I went straightaway to see if there was something to decrease the roundoff error. Off course I didnt get anything of it.

I can send you the whole program so that you can test it and see if it happens the same. My e-mail is guilleh_29@hotmail.com

You know I have a really big doubt about two modules called nrtype and nrutil. I have no idea why do they have to be there. May there be the problem?

About the edition I think it is the first, but I am not sure (I'm doing this at university).

Regards,

Guillermo

...I was told it was better to use Fortran instead of matlab cause it can take it better.
I don't understand. Is the point of the exercise for you to find a solution to a particular set of equations or for you to learn about solving systems of equations that arise from that particular kind of problem?
Matlab is a terrific tool for getting solutions of lots of different things. It is based on some very powerful algorithms developed over the years. With Matlab, you can often get a solution to a system of equations very easily.

(Of course, if the person who told you that it is "better to use Fortran" is an instructor or faculty adviser, further discussion of what to use is irrelevant. Fortran it is.)

However, maybe you can still get some insight from other approaches.

Note that I am using octave (an open-source tool based on the Matlab interface) here; the steps for Matlab would be the same, as far as I know


octave:1> a = [1,2,3; 2,3,4; 1,2,4]
a =

1 2 3
2 3 4
1 2 4

octave:2> b = [14; 20; 17]
b =

14
20
17

octave:3> x = a\b
x =

1
2
3


What happens with your system of equations if you try to get a solution with Matlab? Is there anything unsatisfactory about the "answer" given by Matlab?

Now, if you want some insight into problems with systems of equations like yours, you might want to investigate some of the methods used. Maybe you can find some Matlab documentation that describes the algorithms, or you can investigate the Fortran-based source code for octave. (Lots of luck here---the octave source code is, obviously, not really designed for teaching!).

Or you can read a text like Numerical Recipes that is designed for teaching. Or you can read other texts (or papers or on-line resources) about problems of solving systems of equations.
I am 99% sure it is not a problem related to the matrix.

Now, knowing that your matrix of coefficients is non-singular and knowing the exact solution is one thing. Having a feel for what happens in a real-world set of computations, complete with roundoff error, is something else.



I know how gauss jordan method works, but I have no idea what the gaussj subroutine is exactly doing...I just copied and pasted ...
I think the authors have expressed their intentions that the text (and accompanying code) was designed for teaching purposes, not to be a "cookbook" of subroutines to paste into your applications and blindly call to obtain solutions to specific problems.

For example, using Gauss-Jordan to find an inverse of a matrix and then multiplying the inverse and the right-hand side is not considered a very good strategy for a couple of reasons, although, it can work. If you just want a particular solution to a particular system of 40 equations and 40 unknowns, the run-time efficiency may not be important, but it is generally preferred to use something like LU decomposition (or, perhaps, the computationally-equivalent Gauss Reduction and Back Substitution).

Regardless of the method used to try to obtain a solution, the solution must be considered an approximation, since it may very well be corrupted by roundoff error. People sometimes recommend iterative techniques to try to improve the approximation. Sometimes this helps a lot, and sometimes---not so much. Furthermore, even if plugging in your "solution" into the original equations appears to give something equal to the original right-hand side, the true nature of "the solution" may still be masked by roundoff errors.

Techniques such as SVD (Singular Value Decomposition) may give you a feel for how badly-conditioned your problem is (by showing you the "condition number" of the matrix of coefficients). SVD can (sometimes) lead to a way to minimize the effects of roundoff errors.

Bottom line: Blindly plugging your equation coefficients into some canned routine may give you some numbers, but that approach, by itself, can't give you any insight into the value of the solution or the nature of the problem. At least that's my experience. (Like all of my ramblings, however, that is just my opinion.)


You know I have a really big doubt about two modules called nrtype and nrutil. I have no idea why do they have to be there. May there be the problem?They are there because they define the data types and function/subroutine interfaces that allow routines like gaussj to use other features and functions of the set of routines that make up the Numerical Recipes suite.

About the edition I think it is the first

If it has modules nrtype and nrutil, it is the "Fortran 90" second edition.

Maybe a moderator can move this thread to the appropriate section of the forum.

A final note:
The gaussj subroutine for that version uses single precision arithmetic. Now, simply converting everything to double precision may increase the number of correct significant digits to something acceptable (to you) in the solution for a particular set of equations, or, maybe not. Maybe that's all you need for your project. I can't say. (You still would have to know a little about the program and Fortran, and change the interface specification of the gaussj subroutine so that you could use it with double precision variables. So: how much time have you saved by not having to "program something that was made?" Furthermore, how can you have confidence in a "solution" that is obtained for some system for which you don't know the exact solution? Where's the insight?)

Regards,

Dave

Heys thanks for everything Dave. I decided to apply the LU descomposition instead of Gauss jordan. I found in the book it is requieres a lot less operations.

Regards,

Guillermo