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!

Matrix inverse with LU decomposition in C++

  1. Oct 11, 2014 #1
    1. The problem statement, all variables and given/known data
    the problem is to find the inverse of a 3x3 matrix using LU Decomposition with C++ command, with the numbers designated. in my case, my numbers for the matrix are
    '306
    410
    780'

    Code (Text):
    #include <stdio.h>
    #include <iostream>
    #include <stdlib.h>
    #include <math.h>
    using namespace std;

    int main(void)
    {
        int e, i, j, k, y;
        float A[3][4] = {{3,0,6,1},{4,1,0,2},{7,8,0,3}};
        float x[3][3], c, sum;
        float L[3][4] = {{1,0,0,0},{0,1,0,0},{0,0,1,0}};
            for(j=0; j<=2; j++)
            {
                for(i=0; i<=2; i++)
                {
                    if(i>j)
                    {
                        c=A[i][j]/A[j][j];
                        L[i][j] = c;
                        for(k=0; k<=2; k++)
                        {
                            A[i][k]=A[i][k]-c*A[j][k];
                        }
                    }
                }
            }
                 
            for(y=0; y<=2; y++)
            {
                if(y==0)
                {
                        L[0][3] = 1;
                    L[1][3] = 0;
                    L[2][3] = 0;
                }
                else if (y==1)
                {
                        L[0][3] = 0;
                    L[1][3] = 1;
                    L[2][3] = 0;
                }
                else if (y==2)
                {
                        L[0][3] = 0;
                    L[1][3] = 0;
                    L[2][3] = 1;
                }
                A[0][3]=L[0][3]/L[0][0];
                for(i=1; i<=2; i++)
                {
                    sum=0;
                    for(j=i-1; j>=0; j--)
                    {
                        sum=sum+L[i][j]*A[j][3];
                    }
                    A[i][3]=(L[i][3]-sum)/L[i][i];
                }
                    x[2][y]=A[2][3]/A[2][2];
                    for(i=1; i>=0; i--)
                    {
                        sum=0;
                        for(j=i+1; j<=2; j++)
                        {
                            sum=sum+A[i][j]*x[j][y];
                        }
                        x[i][y]=(A[i][3] - sum)/A[i][i];
                    }
            }
        cout<<"Matriks X"<<endl;
        for(int i=0;i<=2;i++)
        {
            for(int j=0;j<=2;j++)
            {
                cout<< x[i][j]<<"  ";
            }
        cout<<endl;
        }
        return(0);
    }
    2. Relevant equations


    3. The attempt at a solution

    the result from the program will show most numbers correct, except the numbers on X11 and X21

    it's supposed to be
    A11 = A21 = 0
    but it became
    A11 = -3.97364e-08
    A21 = 1.19209e-07

    with loop, i can find out the numbers from each matrix, and all of them are correct

    only the calculations from the X matrix seems wrong
    need help in making A11 and A21 become 0
     
  2. jcsd
  3. Oct 11, 2014 #2

    Mark44

    Staff: Mentor

    Is this the same as the one you posted yesterday? If so, do not repost a problem just because no one has yet responded.

    All of your matrix variables are type float, which can hold only about 6 or 7 decimal places. You will have a lot more precision if you make them double.
     
  4. Oct 11, 2014 #3
    it is a repost, but i think i need to because when i re-read the guidelines, i must use proper English
    which, i think, is not proper on my last post
    so i tried to simplify my words on the new post (here)

    and by using double for all matrix variables, it indeed became more precise, but also became worse
    after changing float -> double
    A11 = 7.40149e-17
    A21 = -2.22045e-16
     
  5. Oct 11, 2014 #4

    Mark44

    Staff: Mentor

    Your English was fine in the other post. If you need to edit a post, and you're no longer able to do it, report the post using the Report button, and a mentor will take care of it.
    Those aren't worse. They're much closer to 0 than the results when you used float variables.

    Floating point arithmetic on computers is inherently imprecise due to roundoff and truncation errors.
     
  6. Oct 11, 2014 #5
    oh, i thought it wasn't fine...
    i see, i'll do so next time

    and also
    although it is closer to 0, it is not 0
    and the actual answer (done manually) is supposed to be 0
    tried changing into integer, but it's even worse than using float or double
    because the solution become numbers without decimals
    is it also because of the rounding and truncation errors?

    and i also notice that there's a minus on A21 and no minus on A11
    a weird change considering the previous one is reversed
     
  7. Oct 11, 2014 #6

    Mark44

    Staff: Mentor

  8. Oct 11, 2014 #7
    Code (Text):
    #include <iomanip>

    cout<< setprecision(3) << fixed << x[i][j]<<"  ";
    (kind of different, but it works the same
    i set it to only spell out 3 digits behind . )

    using that, the result now becomes

    A11 = -0.000
    A21 = 0.000

    well, although it's still weird for 0 to have (-) signs
    but at least better for it shows the (seemingly) same amount as the actual answer

    thanks
     
  9. Oct 11, 2014 #8

    Mark44

    Staff: Mentor

    The setprecision function is just truncating (or possibly rounding - I don't know) the fractional parts of numbers like 0.000000587 and -0.000342, and displaying only the first three decimal places.
     
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: Matrix inverse with LU decomposition in C++
  1. Inverse of a matrix (Replies: 4)

  2. Lu Decomposition in C (Replies: 1)

Loading...