Numerical Method ( Lab Code.c)

Basic Gauss Elimination.cpp

#include<bits/stdc++.h>
#define eps 0.000001

using namespace std;

double a[50][50], b[50], x[50];

void GAUSS(int n, double a[50][50], double b[50], double x[50], int *status)
{
    double pivot, factor, sum;
    /*----------Elimination Starts-------------*/
    for(int k=1; k<=n-1; k++)
    {
        pivot = a[k][k];
        if(pivot<eps)
        {
            *status = 0;
            return;
        }
        *status = 1;
        for(int i=k+1; i<=n; i++)
        {
            factor = a[i][k]/ pivot;
            for(int j=k+1; j<=n; j++)
            {
                a[i][j] = a[i][j] - (factor*a[k][j]);
            }
            b[i] = b[i] - (factor * b[k]);
        }
    }

    /*----------Back Substitution-----------*/

    x[n] = b[n] / a[n][n];
    for(int k=n-1; k>=1; k--)
    {
        sum = 0.0;
        for(int j=k+1; j<=n; j++)
        {
            sum = sum + (a[k][j] * x[j]);
        }
        x[k] = (b[k]-sum) / a[k][k];
    }
    return;

}

int main()
{
    int status, n;

    printf("Number of Equations:\n");
    scanf("%d", &n);

    printf("Co-efficients of N unknowns:\n");
    for(int i=1; i<=n; i++)
        for(int j=1; j<=n; j++)
           scanf("%lf", &a[i][j]);

    printf("Constants of N unknowns:\n");
    for(int i=1; i<=n; i++)
           scanf("%lf", &b[i]);

    GAUSS(n, a, b, x, &status);

    if(status!=0)
    {
        printf("Values of Unknowns are as follows :\n");
        for(int i=1; i<=n; i++)
            printf("%10.6f\n", x[i]);
    }
    else
        printf("Singular Matrix\n");

}

Bisection.cpp

#include<bits/stdc++.h>
#define f(n) (n*n) - 4*n - 10
#define eps 0.00001

using namespace std;

void bisect(double *a, double *b, double *root, int *s, int *cnt)
{
    double x0, x1, x2, f0, f1, f2;
    x1 = *a;
    x2 = *b;
    f1 = f(x1);
    f2 = f(x2);
    if(f1*f2>0)
    {
        *s = 0;
        return;
    }
    else
    {
        *cnt = 0;
        *s = 1;
        while(*cnt<=100)
        {
            x0 = (x1 + x2)/2.0;
            f0 = f(x0);
            if(f0==0.00)
            {
                *root = x0;
                return;
            }
            if((f0 * f1)<0)
            {
                x2 = x0;
            }
            else
            {
                x1 = x0;
                f1 = f0;
            }
            *cnt = *cnt + 1;
             if(abs((x2 - x1)/x2)<eps) break;
        }
        *root = x0;
    }
}

int main()
{
    int cnt = 0 , s;
    double a, b, root;

    printf("Enter the Values of A & B: \n");
    scanf("%lf %lf", &a, &b);

    bisect(&a, &b, &root, &s, &cnt);

    if(s==0)
    {
        printf("Root not found\n");
    }
    else
    {
        printf("Root = %0.16f\n", root);
        printf("F(root) = %f\n", f(root));
        printf("Total Iterations : %d\n", cnt);
    }
}

False Position.cpp

#include<bits/stdc++.h>
#define f(n) (n*n) + n - 2
#define eps 0.00001

using namespace std;

void falsepos(double *a, double *b, double *root, int *s, int *cnt)
{
    double x0, x1, x2, f0, f1, f2;
    x1 = *a;
    x2 = *b;
    f1 = f(x1);
    f2 = f(x2);
    if(f1*f2>0)
    {
        *s = 0;
        return;
    }
    else
    {
        printf("    n        x1            x2\n");
        *cnt = 0;
        *s = 1;
        while(*cnt<=100)
        {
            x0 = x1 - ((f1*(x2-x1))/(f2-f1));
            f0 = f(x0);
            if((f1* f2)<0)
            {
                x2 = x0;
                f2 = f0;
            }
            else
            {
                x1 = x0;
                f1 = f0;
            }
            printf("%5d %10.10lf %15.10lf\n", *cnt, x1, x2);
            *cnt = *cnt + 1;
           if(fabs((x2 - x1)/x2)<eps) break;
        }
        *root = x0;
    }
}

int main()
{
    int cnt = 0 , s;
    double a, b, root;

    printf("Enter the Values of A & B: \n");
    scanf("%lf %lf", &a, &b);

    falsepos(&a, &b, &root, &s, &cnt);

    if(s==0)
    {
        printf("Root not found\n");
    }
    else
    {
        printf("\nRoot = %.10f\n", root);
        printf("F(root) = %.10f\n", f(root));
        printf("Total Iterations : %d\n", cnt);
    }
}

Gauss Elimination With Partial Pivoting.cpp

#include<bits/stdc++.h>
#define eps 0.00001

using namespace std;

double a[50][50], b[50], x[50];

void pivot(int n, double a[50][50], double b[50], int k)
{
    double large, temp;
    int p = k;
    large = fabs(a[k][k]);
    for(int i=k+1; i<=n; i++)
    {
        if(fabs(a[i][k])>large)
        {
            large = fabs(a[i][k]);
            p = i;
        }
    }
    if(p!=k)
    {
        for(int j=k; j<=n; j++)
        {
            temp = a[p][j];
            a[p][j] = a[k][j];
            a[k][j] = temp;
        }

        temp = b[p];
        b[p] = b[k];
        b[k] = temp;
    }
    return;
}

void bsub(int n, double a[50][50], double b[50], double x[50])
{
    x[n] = b[n] / a[n][n];
    double sum;
    for(int k=n-1; k>=1; k--)
    {
        sum = 0.0;
        for(int j=k+1; j<=n; j++)
        {
            sum = sum + (a[k][j] * x[j]);
        }
        x[k] = (b[k]-sum) / a[k][k];
    }
    return;
}

void elim(int n, double a[50][50], double b[50])
{
    double factor;
    for(int k=1; k<=n-1; k++)
    {
        pivot(n,a,b,k);
        for(int i = k+1; i<=n; i++)
        {
            factor = a[i][k]/a[k][k];
            for(int j=k+1; j<=n; j++)
            {
                a[i][j] = a[i][j] - (factor*a[k][j]);
            }
            b[i] = b[i] - (factor * b[k]);
        }
    }
    return;
}

void GAUSS(int n, double a[50][50], double b[50], double x[50])
{
    elim(n, a, b);
    bsub(n,a,b,x);
    return ;
}

int main()
{
    int status, n;

    printf("Number of Equations:\n");
    scanf("%d", &n);

    printf("Co-efficients of N unknowns:\n");
    for(int i=1; i<=n; i++)
        for(int j=1; j<=n; j++)
           scanf("%lf", &a[i][j]);

    printf("Constants of N unknowns:\n");
    for(int i=1; i<=n; i++)
           scanf("%lf", &b[i]);

    GAUSS(n, a, b, x);

    printf("Values of Unknowns are as follows :\n");
    for(int i=1; i<=n; i++)
        printf("%10.6f\n", x[i]);

}

Gauss Seidel.cpp

#include<stdio.h>
#include<math.h>
#define eps 0.000001

int main() {
    int i,j,m,n,l;
    float x[10],a[10][10],b[10],c;
    printf("\nEnter the value of n : \n");
    scanf("%d",&n);
    printf("\nEnter the number of iterations : \n");
    scanf("%d",&l);
    printf("\nEnter the right hand side constants : \n");
    for(i=0;i<n;i++) {
        scanf("%f",&b[i]);
    }
    printf("\nEnter the coefficients row wise : \n");
    for(i=0;i<n;i++) {
        x[i]=0;
        for(j=0;j<n;j++) {
            scanf("%f",&a[i][j]);
        }
    }
    m=1; int deb ; int key;
    line:
      key = 0;
    for(i=0;i<n;i++) {
        c=b[i];
        for(j=0;j<n;j++) {
            if(i!=j) {
                c=c-a[i][j]*x[j];
            }
        }
        double temp = c/a[i][i];
        if(key==0)
        {
            if(fabs((temp-x[i])/temp)>eps)   key = 1;
        }

        x[i]=c/a[i][i];
    }
    m++;
    if(key==1)deb = m;
    if(m<=l) {
        goto line;
    }
    else {
        printf("\nThe Solution is : \n");
        for(i=0;i<n;i++) {
            printf("\nx(%d) = %f\n",i,x[i]);
        }
        printf("\n\nIterations : %d\n", deb);
    }
}

Jacobi Iteration.cpp

#include<bits/stdc++.h>
#define eps 0.000001
#define MAXIT 1000

using namespace std;


void jacobi(int n, float a[50][50], float b[50], float x[50], int *cnt, int *status)
{
    int key;
    float sum, x0[50];
    for(int i=1; i<=n; i++)
        x0[i] = b[i] / a[i][i];

    *cnt = 1;
    lebel:
    key = 0;

    for(int i=1; i<=n; i++)
    {
        sum = b[i];
        for(int j=1; j<=n; j++)
        {
            if(i==j)
                continue;
            sum = sum - a[i][j] * x0[j];
        }
        x[i] = sum / a[i][i];
        if(key==0)
        {
            if(fabs((x[i]-x0[i])/x[i])>eps)   key = 1;
        }
    }
    if(key==1)
    {
        if(*cnt == MAXIT)
        {
            *status = 2;
            return;
        }
        else
        {
            *status = 1;
            for(int i=1; i<=n; i++)   x0[i] = x[i];
        }
        *cnt = *cnt + 1;
        goto lebel;
    }
    return;
}

int main()
{
    int n, cnt, status;
    float a[50][50], b[50], x[50];
    printf("Matrix Size N: ");
    scanf("%d", &n);
    printf("Input CO-efficients:\n");
    for(int i=1; i<=n; i++)
        for(int j=1; j<=n; j++)
            scanf("%f", &a[i][j]);

    printf("Input Vector B:\n");
    for(int i=1; i<=n; i++)  scanf("%f", &b[i]);

    jacobi(n, a, b, x, &cnt, &status);

    if(status==2)
        printf("\nNo Result in %d iterations\n", MAXIT);
    else
    {
        printf("\nSolution Vector X\n");
        for(int i=1; i<=n; i++)
            printf("%15.6f", x[i]);
        printf("\n\nIterations = %d\n", cnt);
    }
}