I am a newbie at C++. I have been writing some code for the past 2 hours. It's about finding a solution to a tridiagonal system.
#include<iostream>
#include<cmath>
using namespace std;
int main()
{
float mat[100][100];
int n;
cout << "Enter the dimention of the matrix: ";
cin >> n;
cout << "Enter the matrix: ";
for (int i = 0; i < n; i++)
{
for (int j = 0; j < n; j++)
{
cin >> mat[i][j];
}
}
float d[100];
cout << "Enter the values of di's: ";
for (int i = 0; i < n; i++)
{
cin >> d[i];
}
for (int i = 0; i < n; i++)
{
cout << "The value of d" << i + 1 << " is" << d[i];
}
cout << " \n";
for (int i = 0; i < n; i++)
{
for (int j = 0; j < n; j++)
cout << mat[i][j] << " ";
{
cout << endl;
}
}
float a[100], b[100], c[100];
for (int i = 0; i < n + 1; i++)
{
b[i] = mat[i][i];
}
for (int l = 0; l < n; l++)
{
cout << "The value of b" << l + 1 << " is " << b[l] << " " << endl;
}
cout << "\n";
for (int j = 0; j < n; j++)
{
a[j] = mat[j][j - 1];
}
for (int i = 1; i < n - 1; i++)
{
cout << "The value of a" << i + 1 << " is: " << a[i] << " " << endl;
}
cout << "\n";
for (int k = 0; k < n - 1; k++)
{
c[k] = mat[k][k + 1];
}
for (int i = 0; i < n - 1; i++)
{
cout << "The value of c" << i + 1 << " is: " << c[i] << " " << endl;
}
//to find alpha
cout << "\n";
float alpha[100];
alpha[0] = b[0];
for (int i = 1; i < n; i++)
{
alpha[i] = b[i] - ((a[i] * c[i - 1]) / alpha[i - 1]);
}
for (int i = 0; i < n; i++)
{
cout << "The value of alpha" << i + 1 << " is: " << alpha[i] << endl;
}
cout << "\n";
//to find beta
float beta[100];
beta[0] = (d[0] / b[0]);
for (int i = 1; i < n; i++)
{
beta[i] = ((d[i] - a[i] * beta[i - 1]) / alpha[i]);
}
for (int i = 0; i < n; i++)
{
cout << " The value of beta" << i + 1 << " is: " << beta[i] << "\n";
}
//finding the solutions
float x[100];
x[n - 1] = beta[n - 1];
for (int i = n - 2; i > 0; i--)
{
x[i] = beta[i] - ((c[i] * x[i + 1]) / alpha[i]);
}
cout << "the solutions are ";
for (int i = 0; i < n - 1; i++)
{
cout << x[i] << "\t";
}
cout << beta[n - 1] << "\t";
}
Output:
PS C:\Users\jitub\desktop> ./a.exe
Enter the dimention of the matrix: 4
Enter the matrix: 3 -1 0 0 -1 3 -1 0 0 -1 3 -1 0 0 -1 3
Enter the values of di's: -5 10 -15 15
The value of d1 is-5The value of d2 is10The value of d3 is-15The value of d4 is15
3 -1 0 0
-1 3 -1 0
0 -1 3 -1
0 0 -1 3
The value of b1 is 3
The value of b2 is 3
The value of b3 is 3
The value of b4 is 3
The value of a2 is: -1
The value of a3 is: -1
The value of c1 is: -1
The value of c2 is: -1
The value of c3 is: -1
The value of alpha1 is: 3
The value of alpha2 is: 2.66667
The value of alpha3 is: 2.625
The value of alpha4 is: 2.61905
The value of beta1 is: -1.66667
The value of beta2 is: 3.125
The value of beta3 is: -4.52381
The value of beta4 is: 4
the solutions are 0 2 -3 4
Using Thomas's algorithm, I have to find the solution. If I do this by hand, the solution I get is {-1, 2, -3, 4}
. But the output I get is always x1=0
, which is not true.
The algorithm says beta4=x4
, and the other values of xi
's are in the for
loop. I have tried many ways to print that value. Still, I couldn't find any way that works.
Also, the value of ai
's should start from a2
. In the last part of this program, I need to print out the values of xi
's from backward.
I think the error is in your implementation of one of the formulas:
Particularly this one:
beta[i] = ((d[i] - a[i] * beta[i - 1]) / alpha[i]);
The correct expression should be (keep reading for another issue):
beta[i] = d[i] - (a[i] * beta[i - 1] / alpha[i]);
This one is also wrong:
x[i] = beta[i] - ((c[i] * x[i + 1]) / alpha[i]);
It should be:
x[i] = (beta[i] - c[i] * x[i + 1]) / alpha[i];
Note how I rearranged the parentheses.
With that fixed, there is something else you may want to check. I don't think you are using the correct indices. It may be just some particularity of your solution but it is my understanding that the expressions should be:
// using a[i-1] instead of a[i]
alpha[i] = b[i] - (a[i-1] * c[i - 1] / alpha[i - 1]);
// using a[1-1] instead of a[i]
beta[i] = d[i] - (a[i-1] * beta[i - 1]) / alpha[i]);
// just for completeness, didn't change the indices
x[i] = (beta[i] - c[i] * x[i + 1]) / alpha[i];