Gaussian Elimination Direct triangular decomposition, single value decomposition, and their error analysis.

Here we will learn about Gaussian Elimination Direct triangular decomposition, single value decomposition, and their error analysis.

Please write the code in C or C++. Please paste the code so that I can run it. Add comments in code so that I can understand better.

Solve the following systems of linear equations using any method and try to estimate the accuracy of the solution.

\begin{pmatrix} 0.000007143 & 0 & 0 & 0\\ 0.6262 & 0.000007355 & 0 & 0\\ 0.4923 & 0.2123 & 0.000002534 & 0\\ 0.8017 & 0.6123 & 0.7165 & 0.000004133 \end{pmatrix}

\[ \begin{pmatrix} x_1 \\ x_2 \\ x_3 \\ x_4 \end{pmatrix} = \begin{pmatrix} .000009245 \\ 0.3763 \\ 0.6087 \\ 0.4306 \end{pmatrix}\]

f_j = f_j-1  + f_j-2 ,   j > 1,   f₀ = 0,  f₁ = 1

(Fibonacci numbers)                      use n = 5, 10, 20, 50

\[ \begin{pmatrix} f_n & f_n+1 \\ f_n+1 & f_n+2 \end{pmatrix} \begin{pmatrix} x_1 \\ x_2 \end{pmatrix} = \begin{pmatrix} f_n+2 \\ f_n+3 \end{pmatrix}\]

Note that for the first matrix the LU decomposition need not be performed. For the first matrix try to solve the equations with another right-hand side vector given by b_i =  \(\sum_{j-1}^{i} \) a_ij and estimate the accuracy.

 

Summary

If the following conditions are met, the matrix is said to be in the reef:

The leading coefficient, or the first non-zero element in each row, is 1.

Each row’s leading coefficient is located in the column to the right of the previous row’s leading coefficient.

Rows with all zeros are at the bottom of the list, followed by rows containing at least one non-zero element.

 

Explanation

If the following conditions are met, the matrix is said to be in r.r.e.f.

All of the requirements for r.e.f.

Each row’s leading coefficient is the column’s only non-zero element.

The method is primarily concerned with performing a series of operations on the matrix’s rows. We want to bear in mind that we want to turn the matrix into an upper triangular matrix in row echelon form while conducting these procedures. The operations can be as follows:

Two rows are swapped.

Using a non-zero scalar to multiply a row

Adding a multiple of one row to another

The procedure is as follows:

Reduction to row echelon structure after forwarding removal. It can be used to determine whether there are no solutions, unique solutions, or an unlimited number of options

Back substitution: a reduction to the reduced row echelon form is made.

 

Code

#include<bits/stdc++.h>
using namespace std;
#define N 3
void backSub(double mat[N][N+1]);
int forwardElim(double mat[N][N+1]);
void gaussianElimination(double mat[N][N+1])
{
    int singular_flag = forwardElim(mat);
    if (singular_flag != -1)
    {
        printf("Singular Matrix.\n");
        if (mat[singular_flag][N])
            printf("Inconsistent System.");
        else
            printf("May have infinitely many solutions.");
        return;
    }
    backSub(mat);
}
void swap_row(double mat[N][N+1], int i, int j)
{
    for (int k=0; k<=N; k++)
    {
        double temp = mat[i][k];
        mat[i][k] = mat[j][k];
        mat[j][k] = temp;
    }
}
void print(double mat[N][N+1])
{
    for (int i=0; i<N; i++, printf("\n"))
        for (int j=0; j<=N; j++)
            printf("%lf ", mat[i][j]);
    printf("\n");
}
int forwardElim(double mat[N][N+1])
{
    for (int k=0; k<N; k++)
    {
        int i_max = k;
        int v_max = mat[i_max][k];
        for (int i = k+1; i < N; i++)
            if (abs(mat[i][k]) > v_max)
                v_max = mat[i][k], i_max = i;
        if (!mat[k][i_max])
            return k; 
        if (i_max != k)
            swap_row(mat, k, i_max);
        for (int i=k+1; i<N; i++)
        {
            double f = mat[i][k]/mat[k][k];
            for (int j=k+1; j<=N; j++)
                mat[i][j] -= mat[k][j]*f;
            mat[i][k] = 0;
        }

    }
    return -1;
}
void backSub(double mat[N][N+1])
{
    double x[N];
    for (int i = N-1; i >= 0; i--)
    {
        x[i] = mat[i][N];
        for (int j=i+1; j<N; j++)
        {
            x[i] -= mat[i][j]*x[j];
        }
        x[i] = x[i]/mat[i][i];
    }

    printf("\nSolution for the system:\n");
    for (int i=0; i<N; i++)
        printf("%lf\n", x[i]);
}
int main()
{
    double mat[N][N+1] =  {{0.000007143,0,0,0},
                       {0.6262,0.000007355,0,0},
                       {0.4923,0.2123,0.000002534,0}
                                          };

    gaussianElimination(mat);

    return 0;
}

Output

 

Also, read find the Intersections Points.