Essential Numerical Methods with C Examples

Essential Numerical Methods with C Examples

1. Newton-Raphson (Nonlinear Equation)

Algorithm (Very Simple Version)

  1. Read starting guess (x)
  2. Read tolerance (how small the error should be)
  3. Repeat many times:
    1. Calculate function value f(x)
    2. Calculate slope f'(x)
    3. If slope is almost zero → stop (bad)
    4. New x = x − f(x) / f'(x)
    5. If change is very small → stop (found answer)
  4. Show the root

C

#include <stdio.h>
#include <math.h>

double f(double x){ return x*x*x - 2*x - 5; }
double df(double x){ return 3*x*x - 2; }

int main(){
    double x, xnew, eps = 0.000001;
    int i = 0;

    printf("Initial guess: "); scanf("%lf",&x);

    do{
        xnew = x - f(x)/df(x);
        printf("x = %.6f\n", xnew);
        if(fabs(xnew - x) < eps) break;
        x = xnew;
        i++;
    }while(i < 50);

    printf("Root ≈ %.6f\n", xnew);
    return 0;
}

2. Trapezoidal Rule (Short)

Algorithm

  1. Read start (a), end (b), number of strips (n)
  2. h = (b – a) / n
  3. sum = f(a) + f(b)
  4. Add 2 × f at every middle point
  5. Result = h/2 × sum

C

#include <stdio.h>

double f(double x){ return 1/(1+x*x); }

int main(){
    double a,b,h,sum=0;
    int n,i;

    printf("a b n: "); scanf("%lf %lf %d",&a,&b,&n);
    h = (b-a)/n;

    sum = f(a) + f(b);
    for(i=1; i<n; i++)
        sum += 2 * f(a + i*h);

    printf("Integral = %.6f\n", h*sum/2);
    return 0;
}

3. Simpson’s 1/3 Rule (Short)

Algorithm

  1. Read a, b, n (n must be even)
  2. h = (b – a) / n
  3. sum = f(a) + f(b)
  4. Add 4 × f(odd points) + 2 × f(even points)
  5. Result = h/3 × sum

C

#include <stdio.h>

double f(double x){ return 1/(1+x*x); }

int main(){
    double a,b,h,sum=0;
    int n,i;

    printf("a b n(even): "); scanf("%lf %lf %d",&a,&b,&n);
    h = (b-a)/n;

    sum = f(a) + f(b);
    for(i=1; i<n; i++)
        sum += (i%2==1 ? 4:2) * f(a+i*h);

    printf("Simpson 1/3 = %.6f\n", h*sum/3);
    return 0;
}

4. Gauss-Seidel Method (Simple)

Algorithm

  1. Read equations (matrix A and b)
  2. Read starting values x1, x2, …
  3. Repeat many times: For each equation → new xi = (b_i – sum of other terms) / a_ii → use new values immediately
  4. Stop when values almost stop changing

C

#include <stdio.h>
#include <math.h>
#define N 3

int main(){
    double a[N][N], b[N], x[N], eps=0.0001;
    int i,j,iter = 0;

    printf("Enter A matrix then b:\n");
    for(i=0;i<N;i++){
        for(j=0;j<N;j++) scanf("%lf",&a[i][j]);
        scanf("%lf",&b[i]);
    }
    printf("Initial x: "); for(i=0;i<N;i++) scanf("%lf",&x[i]);

    do{
        for(i=0;i<N;i++){
            double sum = b[i];
            for(j=0;j<N;j++)
                if(i!=j) sum -= a[i][j]*x[j];
            x[i] = sum / a[i][i];
        }

        printf("x = [%.4f %.4f %.4f]\n",x[0],x[1],x[2]);
        iter++;
    }while(iter<50);   // simple stop (normally check error)

    return 0;
}

5. Gauss Elimination (Basic, No Pivot)

Algorithm

  1. Read matrix A and b (augmented)
  2. For each column k: for each row below k → make element below pivot = 0 by row operation
  3. Back substitution: start from the last equation

C

#include <stdio.h>
#define N 3

int main(){
    double a[N][N+1];
    int i,j,k;

    printf("Enter augmented matrix:\n");
    for(i=0;i<N;i++)
        for(j=0;j<=N;j++)
            scanf("%lf",&a[i][j]);

    // forward elimination
    for(k=0;k<N-1;k++)
        for(i=k+1;i<N;i++){
            double c = a[i][k]/a[k][k];
            for(j=k;j<=N;j++)
                a[i][j] -= c * a[k][j];
        }

    // back substitution
    double x[N];
    for(i=N-1;i>=0;i--){
        x[i] = a[i][N];
        for(j=i+1;j<N;j++)
            x[i] -= a[i][j]*x[j];
        x[i] /= a[i][i];
    }

    printf("Solution:\n");
    for(i=0;i<N;i++) printf("x%d = %.4f\n",i+1,x[i]);
    return 0;
}

6. Runge-Kutta 4 (RK4) — Short Version

Algorithm

  1. Read x₀, y₀, step h, target x
  2. While current x < target:
    1. k1 = h × f(x, y)
    2. k2 = h × f(x + h/2, y + k1/2)
    3. k3 = h × f(x + h/2, y + k2/2)
    4. k4 = h × f(x + h, y + k3)
    5. y = y + (k1 + 2k2 + 2k3 + k4) / 6
    6. x = x + h
  3. Show final y

C

#include <stdio.h>

double f(double x, double y){ return x + y; }

int main(){
    double x,y,h,xn;
    int n,i;

    printf("x0 y0 h xn: "); scanf("%lf %lf %lf %lf",&x,&y,&h,&xn);
    n = (int)((xn-x)/h + 0.5);

    for(i=0; i<n; i++){
        double k1 = h * f(x, y);
        double k2 = h * f(x+h/2, y+k1/2);
        double k3 = h * f(x+h/2, y+k2/2);
        double k4 = h * f(x+h,   y+k3);

        y = y + (k1 + 2*k2 + 2*k3 + k4)/6;
        x += h;

        printf("x=%.3f  y=%.6f\n",x,y);
    }
    return 0;
}