c++numerical-methodslinear-equation

Can someone explain the last part of this code?


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.


Solution

  • 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];