Numerical method : Solution of ordinary differential equation using RK4 method in C

Algorithm:

1. Start
2. Declare and Initialize necessary variable likes K1, K2, K3, K4 etc.
3. Declare and define the function that returns the functional value.
4. Get the initial value of x, initial value of y, no. of iteration and interval from user.
5. for i = 0 to i < n go to step 6.
6. Do the following calculation: and go to step 7
7. The new values of x and y are: x = x + interval, y = y + K;
8. print the value of x and y;
9. stop

Flowchart:

Source Code:

float func(float x, float y){
    return (y*y-x*x)/(y*y+x*x);
}
int main(){
    float K, K1, K2, K3, K4;
    float x0 , y0, x, y;
    int j, n; float i, H;
    printf("Enter initial value of x: ");
    scanf("%f", &x0);
    printf("Enter initial value of y: ");
    scanf("%f", &y0);
    printf("Enter no iteration: ");
    scanf("%d", &n);
    printf("Enter the interval: ");
    scanf("%f", &H);
    x = x0;
    y = y0;
    for(i = x+H, j = 0; j < n; i += H, j++){
        K1 = H * func(x , y);
        K2 = H * func(x+H/2, y+K1/2);
        K3 = H * func(x+H/2, y+K2/2);
        K4 = H * func(x+H, y+K3);
        K = (K1 + 2*K2 + 2*K3 + K4)/6;
        x = i;
        y = y + K;
        printf("At  x = %.2f, y = %.4f ", x, y);
        printf("\n");
    }
    return 0;
}

Output

Numerical Methods: Solution of non-linear equations by using Bisection method in C

Algorithm:

1. Declare and initialize necessary variables like up_range, mid, low_range etc.
2. Read the range( upper and lower)  from user within which the root of the equation  is to be calculated.
3. If root lies within the range? if yes: go to step 4. if no: go to step 2
4. Calculate the mid value of upper and lower range, mid = (upper+lower)/2
5. Calculate the functional value at mid i.e. func(mid).
6. If func(mid)*func(low_range) is less than zero, then replace upper range by mid else replace lower range by mid
7. Display the no of iteration and root
8. if func(mid) is very small? yes: go to step 9. No: go to step 4
9. Display the value of most closest and accurate root.

Source Code:

#include<stdio.h>
#include<math.h>
//function that returns the functional value
float func(float x){
    return (pow(x,3)+5*pow(x,2)-7);
}
int main(){
    float up_range, low_range, mid;
    int i = 0; //no of iteration
    printf("Enter the range: ");
    scanf("%f%f",&up_range,&low_range);
    while(func(up_range)*func(low_range) > 0){ //repeatadly read until the range has root
        printf("\nThis range doesnot contains any root");
        printf("\nEnter again the range: ");
        scanf("%f%f",&up_range,&low_range);
    }
    do{
        mid = (up_range + low_range) / 2;
        if(func(low_range) * func(mid) < 0){ //if signs of mid and low_range is
            up_range = mid;                  //different, replace up_range by mid
        }else{ //else raplace, low_range by mid
            low_range = mid;
        }
        i++;
        printf("\nAt iteration: %d, root = %f",i,mid);
    }while(fabs(func(mid))> 0.0001);
    printf("\nThe root of the equation is %f", mid);
    return 0;
}

Output:

Numerical Methods: Solution of non-linear equation using Newton Raphson method in C

Source Code:

///solution of non-linear  equation using Newton Raphson Method
#include<stdio.h>
#include<math.h>
float f(float x){
    return (pow(x,3) + 5*pow(x,2) -7);
}
float df(float x){
    return (3*pow(x,2) + 10*x );
}
int main(){
    float x0, x1, x2;
    int count = 0;
    printf("Enter the initial guess: ");
    scanf("%f", &x0);
    while(1){
        x1 = x0 - f(x0)/df(x0);
        count++;
        if(x0==x1){
            break;
        }
        x0 = x1;
    }
    printf("The root after %d iteration is %.3f",count, x0);
    return 0;
}

Numerical Methods: Solution of non-linear equation by using Secant method in C

Source Code:

///solution of non linear equations using secant method
#include<stdio.h>
#include<math.h>
float f(float x){
   return (pow(x,5) - 3*pow(x,3)-1);
}
int main(){
    float a, b, c;
    int count = 0;
    printf("Enter the initial value of a: ");
    scanf("%f", &a);
    printf("Enter the initial value of b: ");
    scanf("%f", &b);
    while(1){
         count++;
            c = (a*f(b)-b*f(a))/(f(b)-f(a));
         if(c==b){
            break;
         }
         if(count>=20){
            break;
         }
         a = b;
         b = c;
    }
    printf("\nThe root after %d iteration is %.7f\n",count, c);
    return 0;
}

Numerical Methods: Interpolation with unequal interval with Lagrange’s method in C

Source Code:

///Interpolation with unequal intervals with Lagrange's Formula
#include <stdio.h>
#include <math.h>
#include <windows.h>
COORD coord = {0, 0};
void gotoxy (int x, int y)
{
        coord.X = x; coord.Y = y; // X and Y coordinates
        SetConsoleCursorPosition(GetStdHandle(STD_OUTPUT_HANDLE), coord);
}
int main(){
    float ax[20], ay[20], xp, yp = 0, dr, nr;
    int no_data, i, j;
    printf("How many data? ");
    scanf("%d", &no_data);
    printf("Enter Data: ");
    for(i =0; i < no_data; i++){
        gotoxy(1,2+i); printf("X%d: ", i+1);
        scanf("%f", &ax[i]);
        gotoxy(10, 2+i);printf("Y%d: ", i+1);
        scanf("%f", &ay[i]);
    }
    printf("\nEnter value of x to find corresponding y ");
    scanf("%f", &xp);
    for(i = 0; i < no_data; i++){
        dr = 1;
        nr = 1;
        for(j = 0; j<no_data; j++){
            if(i!=j){
                nr *= xp - ax[j];
                dr *= ax[i] - ax[j];
            }
        }
        yp += nr/dr*ay[i];
    }
    printf("\nrequired value of y is: %f", yp);
    return 0;
}

Numerical Methods:Fitting the curve of the form y = a + bx using least square method in C

///Curve fitting of a straight line of the form y = a + bx
#include<stdio.h>
#include<math.h>
int main(){
    int ax[20], ay[20], i, n;
    int sum_x = 0, sum_xy = 0, sum_x2 = 0, sum_y = 0 ;
    float a , b;
    printf("Enter no of records: ");
    scanf("%d", &n);
    printf("Enter Data: ");
    printf("\nX       Y\n");
    for(i = 0; i < n; i++){
        scanf("%d%d", &ax[i], &ay[i]);
    }
    for(i = 0; i < n; i++){
        sum_x += ax[i];
        sum_y += ay[i];
        sum_xy += ax[i] * ay[i];
        sum_x2 += pow(ax[i], 2);
    }
    b = (n*sum_xy - sum_x*sum_y)/(n*sum_x2 - pow(sum_x,2));
    a = (sum_y - b*sum_x)/n;
    printf("\na = %.3f    b = %.3f", a, b);
    printf("\nThe Equation is %.3f + %.3fX\n\n", a, b);
    return 0;
}

Numerical Methods: Solution of simultaneous algebraic equations using Gauss Elimination method in C

#include<stdio.h>
int main()
{
    double matrix[10][10],a,b, temp[10];
    int  i, j, k, n;
    printf("Enter the no of variables: ");
    scanf("%d", &n);
    printf("Enter the agumented matrix:\n");
    for(i = 0; i < n ; i++){
        for(j = 0; j < (n+1); j++){
            scanf("%lf", &matrix[i][j]);
        }
    }
    for(i = 0; i < n; i++){
       for(j = 0; j < n; j++){
           if(j>i){
                a = matrix[j][i];
                b = matrix[i][i];
                for(k = 0; k < n+1; k++){
                    matrix[j][k] = matrix[j][k] - (a/b) * matrix[i][k];
                }
            }
       }
    }
    printf("The Upper triangular matrix is: \n");
    for( i = 0; i < n; i++){
        for(j = 0; j < n+1; j++){
            printf("%.2f", matrix[i][j]);
            printf("\t");
        }
        printf("\n");
    }
    printf("\nThe required result is: ");
    for(i = n-1; i>=0; i--){
        b = matrix[i][n];
        for(j = n-1 ; j > i; j--){
            b -= temp[n-j]*matrix[i][j];
        }
        temp[n-i] = b/matrix[i][i];
        printf("\n%c => %.2f",97+i, temp[n-i]);
    }
}

Numerical Methods: Solution of simultaneous algebraic equations using Gauss Jordan method in C

#include<stdio.h>
int main()
{
    double matrix[10][10],a,b;
    int  i, j, k, n;
    printf("Enter the no of variables: ");
    scanf("%d", &n);
    printf("Enter the agumented matrix:\n");
    for(i = 0; i < n ; i++){
        for(j = 0; j < (n+1); j++){
            scanf("%lf", &matrix[i][j]);
        }
    }
    for(i = 0; i < n; i++){
       for(j = 0; j < n; j++){
            if(i != j){
                a = matrix[j][i];
                b = matrix[i][i];
                for(k = 0; k < n+1; k++){
                    matrix[j][k] = matrix[j][k] - (a/b) * matrix[i][k];
                }
            }
        }
    }
    for(i = 0; i < n; i++){
         a = matrix[i][i];
        for(j = 0; j < n+1; j++){
            matrix[i][j] /= a;
        }
    }
    printf("The required solution is: \n\n");
    for(i = 0; i < n ; i++){
        printf("%c => %.2f", i+97, matrix[i][n]);
        printf("\n");
    }
    return 0;
}

Numerical Methods: Determinant of nxn matrix using C

Source Code:

#include<stdio.h>
int main(){
    float  matrix[10][10], ratio, det;
    int i, j, k, n;
    printf("Enter order of matrix: ");
    scanf("%d", &n);
    printf("Enter the matrix: \n");
    for(i = 0; i < n; i++){
        for(j = 0; j < n; j++){
            scanf("%f", &matrix[i][j]);
        }
    }
    /* Conversion of matrix to upper triangular */
    for(i = 0; i < n; i++){
        for(j = 0; j < n; j++){
            if(j>i){
                ratio = matrix[j][i]/matrix[i][i];
                for(k = 0; k < n; k++){
                    matrix[j][k] -= ratio * matrix[i][k];
                }
            }
        }
    }
    det = 1; //storage for determinant
    for(i = 0; i < n; i++)
        det *= matrix[i][i];
    printf("The determinant of matrix is: %.2f\n\n", det);
    return 0;
}

Numerical Methods: Inverse of nxn matrix using C

Source Code:

#include<stdio.h>
int main(){
    float matrix[10][10], ratio,a;
    int i, j, k, n;
    printf("Enter order of matrix: ");
    scanf("%d", &n);
    printf("Enter the matrix: \n");
    for(i = 0; i < n; i++){
        for(j = 0; j < n; j++){
            scanf("%f", &matrix[i][j]);
        }
    }
    for(i = 0; i < n; i++){
        for(j = n; j < 2*n; j++){
            if(i==(j-n))
                matrix[i][j] = 1.0;
            else
                matrix[i][j] = 0.0;
        }
    }
    for(i = 0; i < n; i++){
        for(j = 0; j < n; j++){
            if(i!=j){
                ratio = matrix[j][i]/matrix[i][i];
                for(k = 0; k < 2*n; k++){
                    matrix[j][k] -= ratio * matrix[i][k];
                }
            }
        }
    }
    for(i = 0; i < n; i++){
        a = matrix[i][i];
        for(j = 0; j < 2*n; j++){
            matrix[i][j] /= a;
        }
    }
    printf("The inverse matrix is: \n");
    for(i = 0; i < n; i++){
        for(j = n; j < 2*n; j++){
            printf("%.2f", matrix[i][j]);
            printf("\t");
        }
        printf("\n");
    }
    return 0;
}

Numerical Methods: Parabolic curve fitting using C

Source Code:

#include<stdio.h>
#include<math.h>
int main(){
    float xy[20][20], matrix[3][4], ratio, a;
    float sum_x = 0, sum_y = 0, sum_x2 = 0, sum_x3 = 0, sum_x4 = 0, sum_xy = 0, sum_x2y = 0;
    int i, j , k, n;
    printf("Enter no of data: ");
    scanf("%d", &n);
    printf("Enter the data: \n");
    for(i = 0; i < 2; i++){
        for(j = 0; j < n; j++){
            scanf("%f", &xy[i][j]);
        }
    }
    for(i = 0; i < n; i++){
        sum_x += xy[0][i];
        sum_y += xy[1][i];
        sum_x2 += pow(xy[0][i], 2);
        sum_x3 += pow(xy[0][i], 3);
        sum_x4 += pow(xy[0][i], 4);
        sum_xy += xy[0][i]*xy[1][i];
        sum_x2y += pow(xy[0][i], 2) * xy[1][i];
    }
    matrix[0][0] = n;
    matrix[0][1] = sum_x;
    matrix[0][2] = sum_x2;
    matrix[0][3] = sum_y;
    matrix[1][0] = sum_x;
    matrix[1][1] = sum_x2;
    matrix[1][2] = sum_x3;
    matrix[1][3] = sum_xy;
    matrix[2][0] = sum_x2;
    matrix[2][1] = sum_x3;
    matrix[2][2] = sum_x4;
    matrix[2][3] = sum_x2y;
    for(i = 0; i < 3; i++){
        for(j = 0; j < 3; j++){
            if(i!=j){
                ratio = matrix[j][i]/matrix[i][i];
                for(k = 0; k < 4; k++){
                    matrix[j][k] -= ratio * matrix[i][k];
                }
            }
        }
    }
    for(i = 0; i < 3; i++){
        a = matrix[i][i];
        for(j = 0; j < 4; j++){
            matrix[i][j] /= a;
        }
    }
    for(i = 0; i < 3; i++){
        printf("\n%c => %.2f", 97+i, matrix[i][3]);
    }
}

Numerical Methods: Integration of given function using Trapezoidal rule in C

Source Code:

///integration of given function using Trapezoidal rule
#include<stdio.h>
float y(float x){
    return 1/(1+x*x);
}
int main(){
    float x0,xn,h,s;
    int i,n;
    printf("Enter x0, xn, no. of subintervals: ");
    scanf("%f%f%d",&x0,&xn,&n);
    h = (xn-x0)/n;
    s = y(x0) + y(xn);
    for(i = 1; i < n; i++){
        s += 2*y(x0+i*h);
    }
    printf("Value of integral is %6.4f\n",(h/2)*s);
    return 0;
}

Numerical Methods: Integration of given function using Simpson’s 3/8 rule in C

Source Code:

//integration of given function using Simpson's 3/8 rule
#include<stdio.h>
float y(float x){
    return 1/(1+x*x); //function of which integration is to be calculated
}
int main(){
    float x0,xn,h,s;
    int i,n,j,flag;
    printf("Enter x0, xn, no. of subintervals: ");
    scanf("%f%f%d",&x0,&xn,&n);
    h = (xn-x0)/n;
    s = y(x0)+y(xn);
    for(i = 1; i<=n-1;i++){
        for(j=1;j<=n-1;j++){
            if(i==3*j){
                flag = 1;
                break;
            }
            else
                flag = 0;
        }
        if(flag==0)
            s += 3*y(x0+i*h);
        else
            s += 2*y(x0+i*h);
    }
    printf("Value of integral is %6.4f\n",(3*h/8)*s);
    return 0;
}

Numerical Methods: Integration of given function using Simpson’s 1/3 rule in C

Source Code:

///integration of given function using Simpson's 1/3 rule
#include<stdio.h>
float y(float x){
    return 1/(1+x*x);
}
int main(){
    float x0,xn,h,s;
    int i,n;
    printf("Enter x0, xn, no. of subintervals: ");
    scanf("%f%f%d",&x0,&xn,&n);
    h = (xn - x0)/n;
    s = y(x0)+y(xn)+4*y(x0+h);
    for(i = 3; i<=n-1; i+=2){
        s += 4*y(x0+i*h) + 2*y(x0+(i-1)*h);
    }
    printf("Value of integral is %6.4f\n",(h/3)*s);
    return 0;
}

Numerical Methods: Greatest Eigen value and corresponding Eigen vector using power method in C

Source Code:

#include<stdio.h>
#include<math.h>
#include<stdlib.h>
void matrix_mul(float matrix1[3][3],float matrix2[3][1],float matrix3[3][1]){
    int i,j,k;
    for(i = 0; i < 3; i++){
        for(j = 0; j < 1; j++){
            matrix3[i][j] = 0;
            for(k = 0; k < 3; k++){
                matrix3[i][j] += matrix1[i][k]*matrix2[k][j];
            }
        }
    }
}
float findmax(float array[3][1]){
    int i;
    float  maximum;
    maximum = array[0][0];
    for(i = 0; i < 3; i++){
        if(maximum < array[i][0])
            maximum = array[i][0];
    }
    return maximum;
}
int isEqual(float matrix1[3][1], float matrix2[3][1]){
    if(matrix1[0][0] == matrix2[0][0]
    && matrix1[1][0] == matrix2[1][0]
    && matrix1[2][0] == matrix2[2][0])
        return 1;
    return 0;
}
int main(){
    float matrix1[3][3], matrix2[3][1],result[3][1];
    float eigenValue;
    int i,j,k;
    printf("Enter 3x3 matrix:\n");
    for(i = 0; i < 3; i++){
        for(j = 0; j < 3; j++){
            scanf("%f",&matrix1[i][j]);
        }
    }
    for(i = 0; i < 3; i++){
        for(j = 0; j < 1; j++){
            matrix2[i][j] = 1.0;
        }
    }
    while(1){
        matrix_mul(matrix1,matrix2,result);
        eigenValue = findmax(result);
        for(i = 0; i < 3; i++){
            result[i][0]/=eigenValue;
        }
        if(isEqual(matrix2,result)==1){
            break;
        }
        for(i = 0; i < 3; i++){
            matrix2[i][0] = result[i][0];
        }
    }
    printf("Greatest Eigen Value = %f", eigenValue);
    printf("\nAny one of Eigen Vector: \n");
    for(i = 0; i < 3; i++){
        for(j = 0; j < 1; j++){
            printf("%.2f",result[i][j]);
            printf("\t");
        }
        printf("\n");
    }
    return 0;
}