-
Notifications
You must be signed in to change notification settings - Fork 0
/
LU_decompositionForFind_x.c
115 lines (100 loc) · 2.72 KB
/
LU_decompositionForFind_x.c
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
#include <stdio.h>
#define MAX_SIZE 3
void printVector(double vector[MAX_SIZE], int n)
{
for (int i = 0; i < n; i++)
{
printf("%.6lf\n", vector[i]);
}
}
int luDecomposition(double A[MAX_SIZE][MAX_SIZE], double L[MAX_SIZE][MAX_SIZE], double U[MAX_SIZE][MAX_SIZE], int n)
{
for (int i = 0; i < n; i++)
{
for (int k = i; k < n; k++)
{
double sum = 0.0;
for (int j = 0; j < i; j++)
{
sum += L[i][j] * U[j][k];
}
U[i][k] = A[i][k] - sum;
}
for (int k = i; k < n; k++)
{
if (i == k)
{
L[i][i] = 1.0;
}
else
{
double sum = 0.0;
for (int j = 0; j < i; j++)
{
sum += L[k][j] * U[j][i];
}
L[k][i] = (A[k][i] - sum) / U[i][i];
}
}
}
return 1; // Success
}
void forwardSubstitution(double L[MAX_SIZE][MAX_SIZE], double b[MAX_SIZE], double y[MAX_SIZE], int n)
{
for (int i = 0; i < n; i++)
{
y[i] = b[i];
for (int j = 0; j < i; j++)
{
y[i] -= L[i][j] * y[j];
}
y[i] /= L[i][i];
}
}
void backwardSubstitution(double U[MAX_SIZE][MAX_SIZE], double y[MAX_SIZE], double x[MAX_SIZE], int n)
{
for (int i = n - 1; i >= 0; i--)
{
x[i] = y[i];
for (int j = i + 1; j < n; j++)
{
x[i] -= U[i][j] * x[j];
}
x[i] /= U[i][i];
}
}
int main()
{
int n = 3; // Size of the system of equations
double A[MAX_SIZE][MAX_SIZE]; // Coefficient matrix
double b[MAX_SIZE]; // Right-hand side vector
printf("Enter the coefficient matrix A (3x3):\n");
for (int i = 0; i < n; i++)
{
for (int j = 0; j < n; j++)
{
scanf("%lf", &A[i][j]);
}
}
printf("Enter the right-hand side vector b (3x1):\n");
for (int i = 0; i < n; i++)
{
scanf("%lf", &b[i]);
}
double L[MAX_SIZE][MAX_SIZE] = {0}; // Initialize L matrix
double U[MAX_SIZE][MAX_SIZE] = {0}; // Initialize U matrix
double y[MAX_SIZE] = {0}; // Initialize y vector
double x[MAX_SIZE] = {0}; // Initialize x vector
if (luDecomposition(A, L, U, n))
{
forwardSubstitution(L, b, y, n);
backwardSubstitution(U, y, x, n);
printf("Solution (x) Vector:\n");
printVector(x, n);
}
else
{
printf("LU decomposition failed. The system may not have a unique solution.\n");
}
return 0;
}