Numerical Methods
Basic C Programs with Algorithms
Purbanchal University — B.E. (Civil), Year III / Semester II
Basic class-12 style version. No #define macros are used. Each program uses only simple things: a normal
float function for the formula, plain for / while loops, and printf / scanf. Every topic has a
step-by-step Algorithm first, then the C Program, then a Sample Output.
# Program Chapter
1 Bisection Method 2.2
2 Newton-Raphson Method 2.3
3 Secant Method 2.4
4 Fixed Point Iteration Method 2.5
5 Horner's Rule (Polynomial) 2.6
6 Lagrange Interpolation 3.1.2
7 Linear Regression (Least Squares) 3.2.1
8 Trapezoidal Rule 4.2.1
9 Simpson's 1/3 Rule 4.2.2
10 Simpson's 3/8 Rule 4.2.3
11 Basic Gauss Elimination 5.1.1
12 Matrix Inverse (Gauss-Jordan) 5.1.4
13 Euler's Method 6.1
14 Heun's Method 6.2
15 Runge-Kutta (RK4) Method 6.3
Numerical Methods - Basic C Programs Page 1
1. Bisection Method
Chapter 2.2
Algorithm
Step 1: Start.
Step 2: Read interval values a and b.
Step 3: Find midpoint c = (a + b) / 2.
Step 4: If f(a)*f(c) < 0, set b = c, else set a = c.
Step 5: Repeat Steps 3-4 until (b - a) is very small.
Step 6: Print c as the root.
Step 7: Stop.
C Program
#include <stdio.h>
#include <math.h>
// our equation f(x) = x^3 - x - 1
float f(float x)
{
return x*x*x - x - 1;
}
int main()
{
float a, b, c;
printf("Enter a and b: ");
scanf("%f %f", &a, &b);
// keep going until interval is very small
while (fabs(b - a) > 0.0001)
{
c = (a + b) / 2; // midpoint
if (f(a) * f(c) < 0) // root in left part
b = c;
else // root in right part
a = c;
}
printf("Root = %.4f", c);
return 0;
}
Sample Output
Enter a and b: 1 2
Root = 1.3247
Numerical Methods - Basic C Programs Page 2
2. Newton-Raphson Method
Chapter 2.3
Algorithm
Step 1: Start.
Step 2: Read initial guess x.
Step 3: Compute new value x1 = x - f(x) / f'(x).
Step 4: If |x1 - x| is very small, go to Step 6.
Step 5: Set x = x1 and repeat Step 3.
Step 6: Print x1 as the root.
Step 7: Stop.
C Program
#include <stdio.h>
#include <math.h>
// function f(x) = x^3 - x - 1
float f(float x)
{
return x*x*x - x - 1;
}
// derivative f'(x) = 3x^2 - 1
float df(float x)
{
return 3*x*x - 1;
}
int main()
{
float x, x1;
printf("Enter guess: ");
scanf("%f", &x);
do
{
x1 = x - f(x) / df(x); // Newton-Raphson formula
x = x1; // move to new value
}
while (fabs(f(x1)) > 0.0001);
printf("Root = %.4f", x1);
return 0;
}
Sample Output
Enter guess: 1.5
Root = 1.3247
Numerical Methods - Basic C Programs Page 3
3. Secant Method
Chapter 2.4
Algorithm
Step 1: Start.
Step 2: Read two guesses x0 and x1.
Step 3: Compute x2 = x1 - f(x1)*(x1 - x0) / (f(x1) - f(x0)).
Step 4: If |x2 - x1| is very small, go to Step 6.
Step 5: Set x0 = x1, x1 = x2 and repeat Step 3.
Step 6: Print x2 as the root.
Step 7: Stop.
C Program
#include <stdio.h>
#include <math.h>
// function f(x) = x^3 - x - 1
float f(float x)
{
return x*x*x - x - 1;
}
int main()
{
float x0, x1, x2;
printf("Enter x0 and x1: ");
scanf("%f %f", &x0, &x1);
do
{
// secant formula (no derivative needed)
x2 = x1 - f(x1) * (x1 - x0) / (f(x1) - f(x0));
x0 = x1; // shift values
x1 = x2;
}
while (fabs(f(x2)) > 0.0001);
printf("Root = %.4f", x2);
return 0;
}
Sample Output
Enter x0 and x1: 1 2
Root = 1.3247
Numerical Methods - Basic C Programs Page 4
4. Fixed Point Iteration Method
Chapter 2.5
Algorithm
Step 1: Start.
Step 2: Write the equation as x = g(x).
Step 3: Read initial guess x.
Step 4: Compute x1 = g(x).
Step 5: If |x1 - x| is very small, go to Step 7.
Step 6: Set x = x1 and repeat Step 4.
Step 7: Print x1 as the answer.
Step 8: Stop.
C Program
#include <stdio.h>
#include <math.h>
// to find square root of 0.75 we use g(x) = 0.5*(x + 0.75/x)
float g(float x)
{
return 0.5 * (x + 0.75 / x);
}
int main()
{
float x, x1;
printf("Enter guess: ");
scanf("%f", &x);
do
{
x1 = g(x); // x = g(x)
x = x1;
}
while (fabs(g(x1) - x1) > 0.0001);
printf("Answer = %.4f", x1);
return 0;
}
Sample Output
Enter guess: 1
Answer = 0.8660
Numerical Methods - Basic C Programs Page 5
5. Evaluation of Polynomial (Horner's Rule)
Chapter 2.6
Algorithm
Step 1: Start.
Step 2: Read degree n and the (n+1) coefficients.
Step 3: Read the value of x.
Step 4: Set result = first coefficient.
Step 5: For each next coefficient: result = result*x + coefficient.
Step 6: Print result.
Step 7: Stop.
C Program
#include <stdio.h>
int main()
{
int n, i;
float a[20], x, result;
printf("Enter degree: ");
scanf("%d", &n);
printf("Enter coefficients (high power first): ");
for (i = 0; i <= n; i++)
scanf("%f", &a[i]);
printf("Enter x: ");
scanf("%f", &x);
result = a[0]; // start with first coefficient
for (i = 1; i <= n; i++)
result = result * x + a[i]; // Horner step
printf("Value = %.4f", result);
return 0;
}
Sample Output
Polynomial 2x^3 + 3x^2 - 6x + 5 at x = 2
Enter degree: 3
Enter coefficients: 2 3 -6 5
Enter x: 2
Value = 21.0000
Numerical Methods - Basic C Programs Page 6
6. Lagrange Interpolation
Chapter 3.1.2
Algorithm
Step 1: Start.
Step 2: Read number of points n and the x[ ], y[ ] values.
Step 3: Read the value xp where y is needed.
Step 4: Set yp = 0.
Step 5: For each point i, make term = y[i] and multiply it by (xp - x[j]) / (x[i] - x[j]) for all j not equal to i.
Step 6: Add each term to yp.
Step 7: Print yp.
Step 8: Stop.
C Program
#include <stdio.h>
int main()
{
int n, i, j;
float x[20], y[20], xp, yp = 0, term;
printf("Enter number of points: ");
scanf("%d", &n);
printf("Enter x and y values: ");
for (i = 0; i < n; i++)
scanf("%f %f", &x[i], &y[i]);
printf("Enter xp: ");
scanf("%f", &xp);
for (i = 0; i < n; i++)
{
term = y[i]; // start each term with y[i]
for (j = 0; j < n; j++)
{
if (j != i)
term = term * (xp - x[j]) / (x[i] - x[j]);
}
yp = yp + term; // add term to answer
}
printf("Value = %.4f", yp);
return 0;
}
Sample Output
Points (5,12) (6,13) (9,14) (11,16), find y at x=10
Enter number of points: 4
Enter x and y values: 5 12 6 13 9 14 11 16
Enter xp: 10
Value = 14.6667
Numerical Methods - Basic C Programs Page 7
7. Linear Regression (Least Squares)
Chapter 3.2.1
Algorithm
Step 1: Start.
Step 2: Read number of points n and all x, y values.
Step 3: Find sums Sx, Sy, Sxy and Sxx.
Step 4: Compute slope b = (n*Sxy - Sx*Sy) / (n*Sxx - Sx*Sx).
Step 5: Compute intercept a = (Sy - b*Sx) / n.
Step 6: Print the line y = a + b*x.
Step 7: Stop.
C Program
#include <stdio.h>
int main()
{
int n, i;
float x, y, Sx = 0, Sy = 0, Sxy = 0, Sxx = 0, a, b;
printf("Enter number of points: ");
scanf("%d", &n);
for (i = 0; i < n; i++)
{
printf("Enter x and y: ");
scanf("%f %f", &x, &y);
Sx = Sx + x; // sum of x
Sy = Sy + y; // sum of y
Sxy = Sxy + x * y; // sum of x*y
Sxx = Sxx + x * x; // sum of x*x
}
b = (n*Sxy - Sx*Sy) / (n*Sxx - Sx*Sx); // slope
a = (Sy - b*Sx) / n; // intercept
printf("Line: y = %.4f + %.4f x", a, b);
return 0;
}
Sample Output
Data x = 1 2 3 4 5 , y = 3 4 5 6 8
Enter number of points: 5
... (enter each pair)
Line: y = 1.6000 + 1.2000 x
Numerical Methods - Basic C Programs Page 8
8. Trapezoidal Rule
Chapter 4.2.1
Algorithm
Step 1: Start.
Step 2: Read limits a, b and number of strips n.
Step 3: Find strip width h = (b - a) / n.
Step 4: Set sum = f(a) + f(b).
Step 5: For each inside point, add 2*f(x) to sum.
Step 6: Compute answer = (h/2) * sum.
Step 7: Print answer.
Step 8: Stop.
C Program
#include <stdio.h>
// function to integrate f(x) = 1 / (1 + x*x)
float f(float x)
{
return 1.0 / (1.0 + x*x);
}
int main()
{
float a, b, h, x, sum;
int n, i;
printf("Enter a, b and n: ");
scanf("%f %f %d", &a, &b, &n);
h = (b - a) / n; // width of each strip
sum = f(a) + f(b); // two end points
for (i = 1; i < n; i++)
{
x = a + i * h; // inside point
sum = sum + 2 * f(x);
}
printf("Integral = %.4f", (h / 2) * sum);
return 0;
}
Sample Output
Integral of 1/(1+x^2) from 0 to 1
Enter a, b and n: 0 1 10
Integral = 0.7850
Numerical Methods - Basic C Programs Page 9
9. Simpson's 1/3 Rule
Chapter 4.2.2
Algorithm
Step 1: Start.
Step 2: Read limits a, b and n (n must be even).
Step 3: Find h = (b - a) / n.
Step 4: Set sum = f(a) + f(b).
Step 5: For inside points: add 4*f(x) if position is odd, else add 2*f(x).
Step 6: Compute answer = (h/3) * sum.
Step 7: Print answer.
Step 8: Stop.
C Program
#include <stdio.h>
// function to integrate f(x) = 1 / (1 + x*x)
float f(float x)
{
return 1.0 / (1.0 + x*x);
}
int main()
{
float a, b, h, x, sum;
int n, i;
printf("Enter a, b and n (even): ");
scanf("%f %f %d", &a, &b, &n);
h = (b - a) / n;
sum = f(a) + f(b);
for (i = 1; i < n; i++)
{
x = a + i * h;
if (i % 2 == 1) // odd point -> weight 4
sum = sum + 4 * f(x);
else // even point -> weight 2
sum = sum + 2 * f(x);
}
printf("Integral = %.4f", (h / 3) * sum);
return 0;
}
Sample Output
Enter a, b and n (even): 0 1 10
Integral = 0.7854
Numerical Methods - Basic C Programs Page 10
10. Simpson's 3/8 Rule
Chapter 4.2.3
Algorithm
Step 1: Start.
Step 2: Read limits a, b and n (n must be a multiple of 3).
Step 3: Find h = (b - a) / n.
Step 4: Set sum = f(a) + f(b).
Step 5: For inside points: add 2*f(x) if position is a multiple of 3, else add 3*f(x).
Step 6: Compute answer = (3h/8) * sum.
Step 7: Print answer.
Step 8: Stop.
C Program
#include <stdio.h>
// function to integrate f(x) = 1 / (1 + x*x)
float f(float x)
{
return 1.0 / (1.0 + x*x);
}
int main()
{
float a, b, h, x, sum;
int n, i;
printf("Enter a, b and n (multiple of 3): ");
scanf("%f %f %d", &a, &b, &n);
h = (b - a) / n;
sum = f(a) + f(b);
for (i = 1; i < n; i++)
{
x = a + i * h;
if (i % 3 == 0) // every 3rd point -> weight 2
sum = sum + 2 * f(x);
else // others -> weight 3
sum = sum + 3 * f(x);
}
printf("Integral = %.4f", (3 * h / 8) * sum);
return 0;
}
Sample Output
Enter a, b and n (multiple of 3): 0 1 9
Integral = 0.7854
Numerical Methods - Basic C Programs Page 11
11. Basic Gauss Elimination
Chapter 5.1.1
Algorithm
Step 1: Start.
Step 2: Read number of equations n and the augmented matrix.
Step 3: Forward elimination: make all numbers below the diagonal zero using factor = a[i][k] / a[k][k].
Step 4: Back substitution: find the last unknown first, then go upward.
Step 5: Print all unknowns x1, x2, ...
Step 6: Stop.
C Program
#include <stdio.h>
int main()
{
int n, i, j, k;
float a[10][11], x[10], factor;
printf("Enter number of equations: ");
scanf("%d", &n);
printf("Enter augmented matrix: ");
for (i = 0; i < n; i++)
for (j = 0; j <= n; j++)
scanf("%f", &a[i][j]);
// forward elimination
for (k = 0; k < n - 1; k++)
for (i = k + 1; i < n; i++)
{
factor = a[i][k] / a[k][k];
for (j = k; j <= n; j++)
a[i][j] = a[i][j] - factor * a[k][j];
}
// back substitution
for (i = n - 1; i >= 0; i--)
{
x[i] = a[i][n];
for (j = i + 1; j < n; j++)
x[i] = x[i] - a[i][j] * x[j];
x[i] = x[i] / a[i][i];
}
for (i = 0; i < n; i++)
printf("x%d = %.4f\n", i + 1, x[i]);
return 0;
}
Sample Output
2x + y - z = 8 ; -3x - y + 2z = -11 ; -2x + y + 2z = -3
Enter number of equations: 3
Enter augmented matrix: 2 1 -1 8 -3 -1 2 -11 -2 1 2 -3
x1 = 2.0000
x2 = 3.0000
x3 = -1.0000
Numerical Methods - Basic C Programs Page 12
12. Matrix Inverse (Gauss-Jordan)
Chapter 5.1.4
Algorithm
Step 1: Start.
Step 2: Read size n and the matrix.
Step 3: Place an identity matrix on the right side (make it [ A | I ]).
Step 4: For each row: divide the row by its diagonal value to make it 1.
Step 5: Make all other values in that column zero.
Step 6: The right side now becomes the inverse. Print it.
Step 7: Stop.
C Program
#include <stdio.h>
int main()
{
int n, i, j, k;
float a[10][20], pivot, factor;
printf("Enter size: ");
scanf("%d", &n);
printf("Enter matrix: ");
for (i = 0; i < n; i++)
for (j = 0; j < n; j++)
scanf("%f", &a[i][j]);
// make identity matrix on right side
for (i = 0; i < n; i++)
for (j = n; j < 2*n; j++)
{
if (j - n == i)
a[i][j] = 1;
else
a[i][j] = 0;
}
// Gauss-Jordan
for (k = 0; k < n; k++)
{
pivot = a[k][k];
for (j = 0; j < 2*n; j++)
a[k][j] = a[k][j] / pivot; // make diagonal 1
for (i = 0; i < n; i++)
{
if (i != k)
{
factor = a[i][k];
for (j = 0; j < 2*n; j++)
a[i][j] = a[i][j] - factor * a[k][j];
}
}
}
printf("Inverse:\n");
for (i = 0; i < n; i++)
{
for (j = n; j < 2*n; j++)
printf("%.4f ", a[i][j]);
printf("\n");
}
return 0;
}
Numerical Methods - Basic C Programs Page 13
Sample Output
Matrix [2 1 ; 1 1]
Enter size: 2
Enter matrix: 2 1 1 1
Inverse:
1.0000 -1.0000
-1.0000 2.0000
Numerical Methods - Basic C Programs Page 14
13. Euler's Method
Chapter 6.1
Algorithm
Step 1: Start.
Step 2: Read x0, y0, the target xn and step size h.
Step 3: While x < xn, compute y = y + h*f(x, y).
Step 4: Increase x by h each time.
Step 5: Print the final y.
Step 6: Stop.
C Program
#include <stdio.h>
// the equation dy/dx = x + y
float f(float x, float y)
{
return x + y;
}
int main()
{
float x, y, xn, h;
printf("Enter x0, y0, xn, h: ");
scanf("%f %f %f %f", &x, &y, &xn, &h);
while (x < xn)
{
y = y + h * f(x, y); // Euler formula
x = x + h;
}
printf("y(%.2f) = %.4f", xn, y);
return 0;
}
Sample Output
dy/dx = x + y , y(0)=1 , find y(0.2) , h=0.1
Enter x0, y0, xn, h: 0 1 0.2 0.1
y(0.20) = 1.2420
Numerical Methods - Basic C Programs Page 15
14. Heun's Method
Chapter 6.2
Algorithm
Step 1: Start.
Step 2: Read x0, y0, target xn and step size h.
Step 3: Predict yp = y + h*f(x, y).
Step 4: Correct y = y + (h/2)*( f(x, y) + f(x+h, yp) ).
Step 5: Increase x by h. Repeat until x reaches xn.
Step 6: Print the final y.
Step 7: Stop.
C Program
#include <stdio.h>
// the equation dy/dx = x + y
float f(float x, float y)
{
return x + y;
}
int main()
{
float x, y, xn, h, yp;
printf("Enter x0, y0, xn, h: ");
scanf("%f %f %f %f", &x, &y, &xn, &h);
while (x < xn)
{
yp = y + h * f(x, y); // predictor
y = y + (h/2) * (f(x, y) + f(x+h, yp)); // corrector
x = x + h;
}
printf("y(%.2f) = %.4f", xn, y);
return 0;
}
Sample Output
dy/dx = x + y , y(0)=1 , find y(0.2) , h=0.1
Enter x0, y0, xn, h: 0 1 0.2 0.1
y(0.20) = 1.2468
Numerical Methods - Basic C Programs Page 16
15. Runge-Kutta 4th Order (RK4) Method
Chapter 6.3
Algorithm
Step 1: Start.
Step 2: Read x0, y0, target xn and step size h.
Step 3: Find four slopes: k1 = h*f(x, y), k2 = h*f(x+h/2, y+k1/2), k3 = h*f(x+h/2, y+k2/2), k4 = h*f(x+h, y+k3).
Step 4: Update y = y + (k1 + 2k2 + 2k3 + k4)/6.
Step 5: Increase x by h. Repeat until x reaches xn.
Step 6: Print the final y.
Step 7: Stop.
C Program
#include <stdio.h>
// the equation dy/dx = x + y
float f(float x, float y)
{
return x + y;
}
int main()
{
float x, y, xn, h, k1, k2, k3, k4;
printf("Enter x0, y0, xn, h: ");
scanf("%f %f %f %f", &x, &y, &xn, &h);
while (x < xn)
{
k1 = h * f(x, y);
k2 = h * f(x + h/2, y + k1/2);
k3 = h * f(x + h/2, y + k2/2);
k4 = h * f(x + h, y + k3);
y = y + (k1 + 2*k2 + 2*k3 + k4) / 6;
x = x + h;
}
printf("y(%.2f) = %.4f", xn, y);
return 0;
}
Sample Output
dy/dx = x + y , y(0)=1 , find y(0.2) , h=0.1
Enter x0, y0, xn, h: 0 1 0.2 0.1
y(0.20) = 1.2428
Numerical Methods - Basic C Programs Page 17