How could I improve this 3 x 3 matrix linear equation solver?

May 3, 2020 at 3:10pm
Hello,

I have constructed code to solve values for linear equations in matrix form. I have got it working and it calculates correctly compared to my own calculated examples.

What I am really seeking is how I could improve the coding to generate the values for the co factor matrix and calculating determinants for 2x2 and 3x3 ones. Perhaps something involving a for loop or something?

Or maybe there's another way I could be efficient.

By the way, this isn't an assignment. I'm contemplating a career change so I am working through some textbooks.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
// 5.11_Exercises_5.9.cpp : This file contains the 'main' function. Program execution begins and ends there.
/*5.9 Write a module for solving the 3 × 3 linear system Au = b where A is nonsingular*/
    /*  1. Make matrix
        2. Make vector
        3. Evaluate determinant
        4. Form cofactor matrix
        5. Transpose cofactor matrix to obtain adjoint
        6. divide throught the adjoint by the determinant
        7. Multiply A^-1 matrix x b column vector = new column vector for x values*/
#include <iostream>
#include <cmath>

using namespace std;

//Function prototypes for formation and deletion of matrix and vector
double** MakeMatrix(int);
double* MakeVector(int);
void EnterMatrixElements(double**, int);
void DeleteMatrix(double**, int);
void DeleteVector(double*);
void MakeDummyVector(double*);
//Function prototypes for display functions for matrix and vector
void DisplayMatrix(double**, int, int);
void DisplayVector(double*, int);

//Function prototypes to obtain inverse of square matrix functions
double ThreebyThreedeterminant(double**);           // Evaluate determinant of matrix
double TwobyTwodeterminant(double, double, double, double);
double** InverseA(double**);
                                        //Using x = A^-1 . b (b = column vector/matrix formed from linear equation numbers)
double* Mutiply(double**, double*, double*, int, int, int);

double** MakeMatrix(int size)
{
    double** matrix;
    matrix = new double* [size];

    for (int i = 0; i < size; i++)
    {
        matrix[i] = new double[size];
    }
    return matrix;
}
void MakeDummyVector(double* vector)
{
    for (int i = 0; i < 3; i++)
    {
        vector[i] = 0;
    }
}
void EnterMatrixElements(double** matrix, int size)
{
    cout << "What data would you like to enter for the matrix?" << endl;
    cout << "The matrix has a size of " << size * size << " elements" << endl;
    double temp;
    for (int i = 0; i < size; i++)
    {
        for (int j = 0; j < size; j++)
        {
            cout << "What is number you would like to enter?" << endl;
            cin >> temp;
            matrix[i][j] = temp;
        }
    }
    DisplayMatrix(matrix, size, size);
}
double* MakeVector(int size)
{
    //Dynamically allocate memory for vector array
    double* vector = new double[size];
    return vector;
}

void EnterVectorElements(double* vector)
{
    cout << "Enter the values for each element in your column vector" << endl;
    //Enter values for vector array elements
    for (int i = 0; i < 3; i++)
    {
        cin >> vector[i];
    }
}

void DeleteMatrix(double** matrix, int size)
{
    for (int i = 0; i < size; i++)
    {
        delete[] matrix[i];
    }
    delete[] matrix;
}
void DeleteVector(double* vector)
{
    delete[] vector;
}
void DisplayMatrix(double** matrix, int rows, int columns)
{
    cout << "{";
    for (int i = 0; i < rows; i++)
    {
        for (int j = 0; j < columns; j++)
        {
            cout << matrix[i][j] << ", ";
        }
        cout << endl;
    }
    cout << "}";
}
void DisplayVector(double* vector, int size)
{
    cout << "The elements of the vector are {";
    for (int i = 0; i < size; i++)
    {
        cout << vector[i] << endl;
    }
    cout << "}";
}
double TwobyTwodeterminant(double a, double b, double c, double d)
{
    return ((a * d) - (b * c));
}
double ThreebyThreedeterminant(double** A)
{
    double determinant =    A[0][0] * ((A[1][1] * A[2][2]) - (A[1][2] * A[2][1]))
                        -   A[0][1] * ((A[1][0] * A[2][2]) - (A[1][2] * A[2][0]))
                        +   A[0][2] * ((A[1][0] * A[2][1]) - (A[1][1] * A[2][0]));

    return determinant;
}
double** InverseA(double** A) 
{
    double** cofactor = MakeMatrix(3);
    double** cofactor_transposed = MakeMatrix(3);
    double** A_minus_1 = MakeMatrix(3);
    double** tempmatrix = MakeMatrix(3);

    //Form co factor matrix by recursion using the 2 x 2 determinant function using the A matrix
    cofactor[0][0] =      TwobyTwodeterminant(A[1][1], A[1][2], A[2][1], A[2][2]);
    cofactor[0][1] = -1 * TwobyTwodeterminant(A[1][0], A[1][2], A[2][0], A[2][2]);
    cofactor[0][2] =      TwobyTwodeterminant(A[1][0], A[1][1], A[2][0], A[2][1]);
    cofactor[1][0] = -1 * TwobyTwodeterminant(A[0][1], A[0][2], A[2][1], A[2][2]);
    cofactor[1][1] =      TwobyTwodeterminant(A[0][0], A[0][2], A[2][0], A[2][2]);
    cofactor[1][2] = -1 * TwobyTwodeterminant(A[0][0], A[0][1], A[2][0], A[2][1]);
    cofactor[2][0] =      TwobyTwodeterminant(A[0][1], A[0][2], A[1][1], A[1][2]);
    cofactor[2][1] = -1 * TwobyTwodeterminant(A[0][0], A[0][2], A[1][0], A[1][2]);
    cofactor[2][2] =      TwobyTwodeterminant(A[0][0], A[0][1], A[1][0], A[1][1]);

    cout << " The co factor matrix is {" << endl;
    DisplayMatrix(cofactor, 3, 3);

    //Transpose the cofactor matrix
    for (int i = 0; i < 3; i++)
    {
        for (int j = 0; j < 3; j++)
        {
            //Swap columns of cofactors with rows in the traspose
            cofactor_transposed[i][j] = cofactor[j][i];
        }
    }

    cout << " The co factor matrix transposed is " << endl;
    DisplayMatrix(cofactor_transposed, 3, 3);

    //Divide each adjoint element by the determinant of A
    for (int i = 0; i < 3; i++)
    {
        for (int j = 0; j < 3; j++)
        {
            tempmatrix[i][j] = cofactor_transposed[i][j] / ThreebyThreedeterminant(A);
            cout << tempmatrix[i][j] << " = " << cofactor_transposed[i][j] << " / " << ThreebyThreedeterminant(A) << endl;;
        }
    }

    //Assign values to inverse matrix
    for (int i = 0; i < 3; i++)
    {
        for (int j = 0; j < 3; j++)
        {
            A_minus_1[i][j] = tempmatrix[i][j];
        }
    }
    DisplayMatrix(A_minus_1, 3, 3);
    return A_minus_1;
}

double* Mutiply(double** matrix, double* vector1, double* vector2, int rows, int columns, int vector_size)
{
    rows, columns = 3;
    
    cout << "Matrix x Vector " << endl;
    for (int i = 0; i < rows; i++)
    {
        for (int j = 0; j < columns; j++)
        {
            vector2[i] += matrix[i][j] * vector1[j];
        }
    }

    DisplayVector(vector2, 3);

    return vector2;
}

int main()
{
    cout << "This script will solve a system of linear equations in matrix form " << endl;
    double** A = MakeMatrix(3);
    EnterMatrixElements(A, 3);

    double* b = MakeVector(3);
    EnterVectorElements(b);
    double* c = MakeVector(3);
    MakeDummyVector(c);

    cout << ThreebyThreedeterminant(A);

    double** A_Minus_1 = MakeMatrix(3);
    A_Minus_1 = InverseA(A);
    Mutiply(A_Minus_1, b, c, 3, 3, 3);

    DeleteMatrix(A, 3);
    DeleteVector(b);
    DeleteVector(c);
    return 0;
}
/*  1. Make matrix
    2. Make vector
    3. Evaluate determinant
    4. Form cofactor matrix
    5. Transpose cofactor matrix to obtain adjoint
    6. divide throught the adjoint by the determinant
    7. Multiply A^-1 matrix x b column vector = new column vector for x values*/
May 3, 2020 at 7:38pm
Well, since no one else has given input, I might as well say, a for...i, for...j loop wouldn't make it more efficient, it would just save a bit of code real estate. As I'm sure you know, the (i)(j)th minor is obtained by taking the det of the matrix without the i'th row and j'th column. I'm sure if I thought hard enough about it, I could turn it into a generalized loop (this would be necessary in the general NxN case). But if you only want to deal with a 2x2 or 3x3 matrix, having the loop "unrolled" as you do is fine.

Maybe Gaussian elimination method would be faster than Cramer's rule for large matrices, but I sincerely doubt it matters for such small sizes.

From just a code maintainability point of view, I would avoid printing things in the functions where you compute things. For example, I would have InverseA simply return the inverse of A. If you want to print the result out, do it outside of the computation function itself. This just helps code re-use/maintainability a bit.
Last edited on May 3, 2020 at 7:46pm
May 6, 2020 at 1:36pm
Hi,

I can see where you are coming from for the generalised case. Incidentally, I had to calculate values for N x N systems of equations! Also, from a code re-use standpoint, I guess that having the function there to calculate and not display avoids conflicts with recursively calling display functions in different programs? Since you'd want to display in a particular way. I just did it so I could check values as I went along.
Topic archived. No new replies allowed.