Matrix inverse with LU decomposition in C++

Click For Summary
SUMMARY

The discussion focuses on finding the inverse of a 3x3 matrix using LU decomposition in C++. The provided code initially uses float data types for matrix variables, resulting in imprecise values for certain elements of the output matrix. After switching to double data types, the precision improved, but the output still did not yield exact zeros for specific elements, highlighting the inherent limitations of floating-point arithmetic. The use of the setprecision function from the iostream library was suggested to format the output more clearly, allowing for better presentation of results.

PREREQUISITES
  • Understanding of LU decomposition for matrix inversion
  • Proficiency in C++ programming, specifically with arrays and loops
  • Familiarity with floating-point arithmetic and its limitations
  • Knowledge of the iostream library and formatting output in C++
NEXT STEPS
  • Explore advanced matrix operations using Eigen library in C++
  • Learn about numerical stability in floating-point computations
  • Investigate the implications of using different data types (float vs. double) in C++
  • Study the implementation of matrix operations in other programming languages, such as Python with NumPy
USEFUL FOR

Students and professionals in computer science, particularly those focusing on numerical methods, linear algebra, and C++ programming for scientific computing.

KWKWII
Messages
4
Reaction score
0

Homework Statement


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:
#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);
}

Homework Equations

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
 
Physics news on Phys.org
Is this the same as the one you posted yesterday? If so, do not repost a problem just because no one has yet responded.

KWKWII said:

Homework Statement


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:
#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);
}

Homework Equations

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
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.
 
Mark44 said:
do not repost a problem just because no one has yet responded

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
 
KWKWII said:
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)
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.
KWKWII said:
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
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.
 
Mark44 said:
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

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
 
Mark44 said:
using the precision function

Code:
#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
 
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.
 

Similar threads

  • · Replies 7 ·
Replies
7
Views
2K
  • · Replies 17 ·
Replies
17
Views
2K
  • · Replies 10 ·
Replies
10
Views
3K
  • · Replies 13 ·
Replies
13
Views
2K
  • · Replies 3 ·
Replies
3
Views
1K
  • · Replies 13 ·
Replies
13
Views
3K
  • · Replies 14 ·
Replies
14
Views
4K
  • · Replies 5 ·
Replies
5
Views
3K
  • · Replies 23 ·
Replies
23
Views
9K
  • · Replies 10 ·
Replies
10
Views
2K