Linear Regression

C code:

#include<stdio.h>
#include<conio.h>
int main()
{
int n,i,j,k;
float a = 0, b = 0 ,x[10] ,y[10], sx = 0 ,sy = 0 ,sxy=0, sx2 = 0;
printf("Enter the number of points \n") ;
scanf("%d",&n);
printf("Enter the value of x and fx \n");
for(i=0;i<n;i++)
{
scanf("%f%f",&x[i],&y[i]);
}
for(i=0;i<n;i++)
{
sx = sx + x[i];
sy = sy + y[i];
sxy = sxy+x[i]*y[i];
sx2 = sx2+x[i] * x[i] ;
}
b=((n*sxy) - (sx*sy))/((n*sx2) - (sx*sx));
a= (sy / n) -(b * sx/n);
printf("Fitted line is: %f + % f x ",a,b);
getch();
}

Python code:

import numpy as np

def linear_regression(X, y):

return np.linalg.lstsq(np.c_[np.ones(len(X)), X], y, rcond=None)[0]

X = np.array([list(map(float, input().split()))

for _ in range(int(input("Enter number of samples: ")))])

y = np.array(list(map(float, input().split())))

print("Linear Regression Coefficients:", linear_regression(X, y))

Polynomial Regression

C code:

#include <stdio.h>

#include <math.h>

#define N 5

#define DEG 2

void solve(double A[DEG+1][DEG+2], int n) {

for (int i = 0; i <= n; i++) {

for (int k = i+1; k <= n; k++) {

double ratio = A[k][i] / A[i][i];

for (int j = i; j <= n+1; j++)

A[k][j] -= ratio * A[i][j];

double coeffs[DEG+1];

for (int i = n; i >= 0; i--) {

coeffs[i] = A[i][n+1];

for (int j = i+1; j <= n; j++)

coeffs[i] -= A[i][j] * coeffs[j];

coeffs[i] /= A[i][i];

printf("Polynomial Coefficients:\n");

for (int i = 0; i <= n; i++) printf("a%d = %.4f\n", i, coeffs[i]);

int main() {

double x[N] = {1, 2, 3, 4, 5};

 double y[N] = {2.2, 2.8, 3.6, 4.5, 5.1};

double X[2*DEG+1] = {0}, B[DEG+1] = {0}, A[DEG+1][DEG+2] = {0};

for (int i = 0; i < N; i++)

for (int j = 0; j <= 2*DEG; j++)

X[j] += pow(x[i], j);

for (int i = 0; i <= DEG; i++) {

for (int j = 0; j <= DEG; j++)

A[i][j] = X[i+j];

for (int k = 0; k < N; k++)

B[i] += pow(x[k], i) * y[k];

A[i][DEG+1] = B[i];

solve(A, DEG);

return 0;

Python code:

import numpy as np

x = np.array([1, 2, 3, 4, 5])
y = np.array([2.2, 2.8, 3.6, 4.5, 5.1])

coeffs = np.polyfit(x, y, deg=2)

print("Polynomial Coefficients:", coeffs)

Exponential Regression

C code:

#include <stdio.h>

#include <math.h>

void exponential_regression(int n, double x[], double y[], double *a, double *b) {

double sum_x = 0, sum_y = 0, sum_xy = 0, sum_x2 = 0;

for (int i = 0; i < n; i++) {

sum_x += x[i];

sum_y += log(y[i]);

sum_xy += x[i] * log(y[i]);

sum_x2 += x[i] * x[i];

double denom = n * sum_x2 - sum_x * sum_x;

*b = (n * sum_xy - sum_x * sum_y) / denom;

*a = exp((sum_y - (*b) * sum_x) / n);

int main() {

double x[] = {1, 2, 3, 4, 5}, y[] = {2.2, 2.8, 3.6, 4.5, 5.1};

int n = sizeof(x) / sizeof(x[0]);

double a, b;

exponential_regression(n, x, y, &a, &b);

printf("y = %.4f * e^(%.4fx)\n", a, b);

return 0;

Python code:

import numpy as np
def exponential_regression(x, y):

log_y = np.log(y)

b, log_a = np.polyfit(x, log_y, 1)

return np.exp(log_a), b

x = np.array([1, 2, 3, 4, 5])

y = np.array([2.2, 2.8, 3.6, 4.5, 5.1])

a, b = exponential_regression(x, y)

print(f"y = {a:.4f} * e^({b:.4f}x)")

Trapezoidal Rule & Composite Trapezoidal Rule

C code:

#include <stdio.h>

#include <math.h>

double f(double x) {

return x * x;

double trapezoidal(double a, double b) {

return (b - a) * (f(a) + f(b)) / 2.0;

double composite_trapezoidal(double a, double b, int n) {

double h = (b - a) / n, sum = f(a) + f(b);

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

return (h / 2) * sum;

int main() {

double a = 0, b = 1;

int n = 4;

printf("Trapezoidal Rule: %lf\n", trapezoidal(a, b));

printf("Composite Trapezoidal Rule: %lf\n", composite_trapezoidal(a, b, n));

return 0;

Python code:
def f(x):

return x**2

def trapezoidal(a, b):

return (b - a) * (f(a) + f(b)) / 2

def composite_trapezoidal(a, b, n):

h = (b - a) / n

return (h / 2) * (f(a) + f(b) + 2 * sum(f(a + i * h) for i in range(1, n)))

a, b, n = 0, 1, 4

print("Trapezoidal Rule:", trapezoidal(a, b))

print("Composite Trapezoidal Rule:", composite_trapezoidal(a, b, n))

Simpson’s 1/3 Rule & Composite Simpson’s 1/3 Rule

C code:

#include <stdio.h>

#include <math.h>

double f(double x) {

return x * x;

double simpson13(double a, double b) {

double h = (b - a) / 2.0;

return (h / 3) * (f(a) + 4 * f(a + h) + f(b));

double compositeSimpson13(double a, double b, int n) {

if (n % 2 != 0) n++;

double h = (b - a) / n, sum = f(a) + f(b);

for (int i = 1; i < n; i++)

sum += (i % 2 ? 4 : 2) * f(a + i * h);

return (h / 3) * sum;

int main() {

double a = 0, b = 1;

int n = 4;

printf("Simpson's 1/3 Rule: %.6lf\n", simpson13(a, b));

printf("Composite Simpson's 1/3 Rule: %.6lf\n", compositeSimpson13(a, b, n));

return 0;
}

Python code:

import numpy as np

def f(x):

return x ** 2

def simpson13(a, b):

h = (b - a) / 2.0

return (h / 3) * (f(a) + 4 * f(a + h) + f(b))

def composite_simpson13(a, b, n):

if n % 2 != 0:

n += 1

h = (b - a) / n

x = np.linspace(a, b, n + 1)

return (h / 3) * (f(x[0]) + 4 * sum(f(x[1:n:2])) + 2 * sum(f(x[2:n-1:2])) + f(x[n]))

a, b, n = 0, 1, 4

print("Simpson's 1/3 Rule:", simpson13(a, b))

print("Composite Simpson's 1/3 Rule:", composite_simpson13(a, b, n))

Simpson’s 3/8 rule

C code:

#include <stdio.h>

#include <math.h>

double f(double x) {

return pow(x, 3);

double simpsons_38(double a, double b) {

double h = (b - a) / 3.0;

return (3.0 * h / 8.0) * (f(a) + 3*f(a+h) + 3*f(a+2*h) + f(b));

double composite_simpsons_38(double a, double b, int n) {

if (n % 3 != 0) n += (3 - n % 3);

double h = (b - a) / n, sum = f(a) + f(b);

for (int i = 1; i < n; i++)

sum += (i % 3 == 0 ? 2 : 3) * f(a + i * h);

return (3 * h / 8) * sum;

int main() {

double a = 1, b = 2;

int n = 6;

printf("Simpson's 3/8 Rule: %.6f\n", simpsons_38(a, b));

 printf("Composite Simpson's 3/8 Rule: %.6f\n", composite_simpsons_38(a, b, n));

return 0;

Python code:

def f(x):

return x**3

def simpsons_38(a, b):

h = (b - a) / 3

return (3 * h / 8) * (f(a) + 3*f(a+h) + 3*f(a+2*h) + f(b))

def composite_simpsons_38(a, b, n):

if n % 3 != 0:

n += (3 - n % 3)

h = (b - a) / n

sum_ = f(a) + f(b) + sum((3 if i % 3 else 2) * f(a + i * h) for i in range(1, n))

return (3 * h / 8) * sum_

a, b, n = 1, 2, 6

print("Simpson's 3/8 Rule:", simpsons_38(a, b))

print("Composite Simpson's 3/8 Rule:", composite_simpsons_38(a, b, n))

Gauss Elimination Method

C code:

#include <stdio.h>

#include <math.h>

int main() {

int n, i, j, k;

float a[10][10], b[10], x[10], factor;

printf("Enter dimension of system of equations: ");

scanf("%d", &n);

printf("Enter coefficients row-wise:\n");

for (i = 0; i < n; i++)

for (j = 0; j < n; j++)

scanf("%f", &a[i][j]);

printf("Enter RHS vector:\n");

for (i = 0; i < n; i++)

scanf("%f", &b[i]);

for (k = 0; k < n - 1; k++) {

if (fabs(a[k][k]) < 1e-6) {

printf("Method failed\n");

return 1;

for (i = k + 1; i < n; i++) {

factor = a[i][k] / a[k][k];

 for (j = k; j < n; j++)

a[i][j] -= factor * a[k][j];

b[i] -= factor * b[k];

for (i = n - 1; i >= 0; i--) {

x[i] = b[i];

for (j = i + 1; j < n; j++)

x[i] -= a[i][j] * x[j];

x[i] /= a[i][i];

for (i = 0; i < n; i++)

printf("x%d = %f\n", i + 1, x[i]);

return 0;

Python code:

import numpy as np

def gaussian_elimination(A, b):

n = len(b)

A, b = A.astype(float), b.astype(float)

for k in range(n - 1):

if abs(A[k, k]) < 1e-6:

return "Method failed"

for i in range(k + 1, n):

 factor = A[i, k] / A[k, k]

A[i, k:] -= factor * A[k, k:]

b[i] -= factor * b[k]

x = np.zeros(n)

for i in range(n - 1, -1, -1):

x[i] = (b[i] - A[i, i + 1:] @ x[i + 1:]) / A[i, i]

return x

n = int(input("Enter dimension: "))

A = np.array([list(map(float, input().split())) for _ in range(n)])

b = np.array(list(map(float, input().split())))

print("Gaussian Elimination Solution:", gaussian_elimination(A, b))

Gauss Jordan Method

C output:

#include <stdio.h>

void gaussJordan(float matrix[3][4]) {

for (int i = 0; i < 3; i++) {

float diag = matrix[i][i];

for (int j = 0; j < 4; j++) matrix[i][j] /= diag;

for (int k = 0; k < 3; k++) {

if (k != i) {

float factor = matrix[k][i];

for (int j = 0; j < 4; j++) matrix[k][j] -= factor * matrix[i][j];

int main() {

float matrix[3][4] = {

{2, 1, -1, 8},

{-3, -1, 2, -11},

{-2, 1, 2, -3}

};

gaussJordan(matrix);

for (int i = 0; i < 3; i++) {

printf("x%d = %.2f\n", i + 1, matrix[i][3]);

}
 return 0;

Python code:

def gauss_jordan(matrix):

for i in range(len(matrix)):

diag = matrix[i][i]

matrix[i] = [x / diag for x in matrix[i]]

for j in range(len(matrix)):

if i != j:

factor = matrix[j][i]

matrix[j] = [matrix[j][k] - factor * matrix[i][k] for k in range(len(matrix[i]))]

matrix = [

[2, 1, -1, 8],

[-3, -1, 2, -11],

[-2, 1, 2, -3]

gauss_jordan(matrix)

for i in range(len(matrix)):

print(f"x{i+1} = {matrix[i][3]:.2f}")
Matrix Inversion by Gauss Jordan Method

C code:

#include <stdio.h>

#define N 3

void gaussJordan(float a[N][N], float inv[N][N]) {

for (int i = 0; i < N; i++)

for (int j = 0; j < N; j++)

inv[i][j] = (i == j);

for (int i = 0; i < N; i++) {

float pivot = a[i][i];

for (int j = 0; j < N; j++) {

a[i][j] /= pivot;

inv[i][j] /= pivot;

for (int k = 0; k < N; k++) {

if (k != i) {

float factor = a[k][i];

for (int j = 0; j < N; j++) {

a[k][j] -= factor * a[i][j];

inv[k][j] -= factor * inv[i][j];

int main() {
 float a[N][N] = {{2, 1, 1}, {1, 3, 2}, {1, 0, 0}}, inv[N][N];

gaussJordan(a, inv);

for (int i = 0; i < N; i++, printf("\n"))

for (int j = 0; j < N; j++)

printf("%.2f ", inv[i][j]);

return 0;

Python code:

import numpy as np

def gauss_jordan(a):

n = len(a)

a = np.array(a, dtype=float)

inv = np.eye(n)

for i in range(n):

a[i] /= a[i, i]

inv[i] /= a[i, i]

for k in range(n):

if k != i:

factor = a[k, i]

a[k] -= factor * a[i]

inv[k] -= factor * inv[i]

return inv

A = [[2, 1, 1], [1, 3, 2], [1, 0, 0]]

inv = gauss_jordan(A)

print(np.round(inv, 2))
Matrix Factorization LU Method

C code:

#include <stdio.h>

void luDecomposition(int n, float a[][10], float l[][10], float u[][10]) {

for (int i = 0; i < n; i++) {

for (int k = i; k < n; k++) {

u[i][k] = a[i][k];

for (int j = 0; j < i; j++) u[i][k] -= l[i][j] * u[j][k];

for (int k = i; k < n; k++) {

if (i == k) l[i][i] = 1;

else {

l[k][i] = a[k][i];

for (int j = 0; j < i; j++) l[k][i] -= l[k][j] * u[j][i];

l[k][i] /= u[i][i];

int main() {

int n;

float a[10][10], l[10][10] = {0}, u[10][10] = {0};

printf("Enter matrix size: ");

scanf("%d", &n);

printf("Enter matrix elements:\n");

for (int i = 0; i < n; i++)

 for (int j = 0; j < n; j++) scanf("%f", &a[i][j]);

luDecomposition(n, a, l, u);

printf("Lower Matrix:\n");

for (int i = 0; i < n; i++, printf("\n"))

for (int j = 0; j < n; j++) printf("%0.2f ", l[i][j]);

printf("Upper Matrix:\n");

for (int i = 0; i < n; i++, printf("\n"))

for (int j = 0; j < n; j++) printf("%0.2f ", u[i][j]);

return 0;

Python code:

def lu_decomposition(n, a):

l = [[0.0 for _ in range(n)] for _ in range(n)]

u = [[0.0 for _ in range(n)] for _ in range(n)]

for i in range(n):

for k in range(i, n):

u[i][k] = a[i][k]

for j in range(i):

u[i][k] -= l[i][j] * u[j][k]

for k in range(i, n):

if i == k:
 l[i][i] = 1.0

else:

l[k][i] = a[k][i]

for j in range(i):

l[k][i] -= l[k][j] * u[j][i]

l[k][i] /= u[i][i]

return l, u

def main():

n = int(input("Enter matrix size: "))

print("Enter matrix elements:")

a = []

for i in range(n):

row = list(map(float, input().split()))

a.append(row)

l, u = lu_decomposition(n, a)

print("Lower Matrix:")

for row in l:

print(" ".join(f"{val:.2f}" for val in row))

print("Upper Matrix:")

for row in u:
 print(" ".join(f"{val:.2f}" for val in row))

if __name__ == "__main__":

main()
Jacobi Iterative Method

C code:

#include <stdio.h>

#include <math.h>

#include <stdlib.h>

void jacobi(int n, double A[10][10], double b[10], double tol, int max_iter, double x[10]) {

double x_new[10];

int iter;

for (int i = 0; i < n; i++) {

x[i] = 0.0;

for (iter = 0; iter < max_iter; iter++) {

for (int i = 0; i < n; i++) {

double sum = 0.0;

for (int j = 0; j < n; j++) {

if (j != i) {

sum += A[i][j] * x[j];

x_new[i] = (b[i] - sum) / A[i][i];

double norm = 0.0;

for (int i = 0; i < n; i++) {

 if (fabs(x_new[i] - x[i]) > norm) {

norm = fabs(x_new[i] - x[i]);

if (norm < tol) {

for (int i = 0; i < n; i++) {

x[i] = x_new[i];

return;

for (int i = 0; i < n; i++) {

x[i] = x_new[i];

int main() {

int n, max_iter = 100;

double tol = 1e-6;

double A[10][10], b[10], x[10];

printf("Enter dimension: ");

scanf("%d", &n);

printf("Enter matrix A (row-wise):\n");

 for (int i = 0; i < n; i++) {

for (int j = 0; j < n; j++) {

scanf("%lf", &A[i][j]);

printf("Enter vector b:\n");

for (int i = 0; i < n; i++) {

scanf("%lf", &b[i]);

jacobi(n, A, b, tol, max_iter, x);

printf("Jacobi Iterative Solution:\n");

for (int i = 0; i < n; i++) {

printf("x[%d] = %.6f\n", i, x[i]);

return 0;

Python code:

import numpy as np

def jacobi(A, b, tol=1e-6, max_iter=100):

n, x = len(b), np.zeros_like(b, dtype=float)

for _ in range(max_iter):
 x_new = (b - (A @ x - np.diag(A) * x)) / np.diag(A)

if np.linalg.norm(x_new - x, np.inf) < tol:

return x_new

x = x_new

return x

n = int(input("Enter dimension: "))

A = np.array([list(map(float, input().split())) for _ in range(n)])

b = np.array(list(map(float, input().split())))

print("Jacobi Iterative Solution:", jacobi(A, b))

Gauss Seidel Method

C code:

#include <stdio.h>

#include <math.h>

void gaussSeidel(int n, double A[n][n], double b[n], double x[n], double tol, int max_iter) {

for (int iter = 0; iter < max_iter; iter++) {

double max_diff = 0;

for (int i = 0; i < n; i++) {

double sum = b[i];

for (int j = 0; j < n; j++)

if (i != j) sum -= A[i][j] * x[j];

double x_new = sum / A[i][i];

max_diff = fmax(max_diff, fabs(x_new - x[i]));

x[i] = x_new;

if (max_diff < tol) return;

int main() {

int n = 4;

double A[4][4] = {{4, -1, 0, 0}, {-1, 4, -1, 0}, {0, -1, 4, -1}, {0, 0, -1, 3}};

double b[4] = {15, 10, 10, 10}, x[4] = {0, 0, 0, 0};

gaussSeidel(n, A, b, x, 1e-6, 100);

printf("Solution:\n");

for (int i = 0; i < n; i++) printf("x%d = %f\n", i+1, x[i]);

 return 0;

Python code:

import numpy as np

def gauss_seidel(A, b, tol=1e-6, max_iter=100):

n, x = len(A), np.zeros_like(b, dtype=float)

for _ in range(max_iter):

x_new = np.copy(x)

for i in range(n):

x_new[i] = (b[i] - np.dot(A[i, :i], x_new[:i]) - np.dot(A[i, i+1:], x[i+1:])) / A[i, i]

if np.linalg.norm(x_new - x, ord=np.inf) < tol:

return x_new

x = x_new

return x

A = np.array([[4, -1, 0, 0], [-1, 4, -1, 0], [0, -1, 4, -1], [0, 0, -1, 3]], dtype=float)

b = np.array([15, 10, 10, 10], dtype=float)

print("Solution:", gauss_seidel(A, b))

Euler’s Method

C output:

#include <stdio.h>

#include <math.h>

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

double euler(double (*f)(double, double), double x0, double y0, double h, int n) {

for (int i = 0; i < n; i++) {

y0 += h * f(x0, y0);

x0 += h;

return y0;

int main() {

double x0, y0, h;

int n;

printf("Enter initial x, y: ");

scanf("%lf %lf", &x0, &y0);

printf("Enter step size h and number of steps: ");

scanf("%lf %d", &h, &n);

printf("Euler's Method Solution: %lf\n", euler(f, x0, y0, h, n));

return 0;

Python output:

def euler(f, x0, y0, h, n):

for _ in range(n):
y0 += h * f(x0, y0)
x0 += h
return y0

expr = input("Enter dy/dx function (in terms of x and y): ")

f = lambda x, y: eval(expr)

x0, y0 = map(float, input("Enter initial x, y: ").split())

h, n = map(int, input("Enter step size h and number of steps: ").split())

print("Euler's Method Solution:", euler(f, x0, y0, h, n))

Heun’s Method

C code:

#include <stdio.h>

double heun(double (*f)(double, double), double x0, double y0, double h, double xn) {

while (x0 < xn) {

double y_pred = y0 + h * f(x0, y0);

y0 += (h / 2) * (f(x0, y0) + f(x0 + h, y_pred));

x0 += h;

return y0;

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

int main() {

printf("%lf\n", heun(f, 0, 1, 0.1, 0.5));

return 0;

Python code:

import numpy as np

def heun(f, x0, y0, h, xn):

while x0 < xn:
y_pred = y0 + h * f(x0, y0)
y0 += (h / 2) * (f(x0, y0) + f(x0 + h, y_pred))
x0 += h
return y0

f = lambda x, y: x + y
print(heun(f, 0, 1, 0.1, 0.5))
Shooting Method

C code:

#include <stdio.h>

#include <math.h>

#define TOL 1e-6

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

double g(double x, double y, double z) { return -y; }

void rk4(double *y, double *z, double x, double h) {

double k1 = h * f(x, *y, *z), l1 = h * g(x, *y, *z);

double k2 = h * f(x + h / 2, *y + k1 / 2, *z + l1 / 2);

double l2 = h * g(x + h / 2, *y + k1 / 2, *z + l1 / 2);

double k3 = h * f(x + h / 2, *y + k2 / 2, *z + l2 / 2);

double l3 = h * g(x + h / 2, *y + k2 / 2, *z + l2 / 2);

double k4 = h * f(x + h, *y + k3, *z + l3);

double l4 = h * g(x + h, *y + k3, *z + l3);

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

*z += (l1 + 2*l2 + 2*l3 + l4) / 6;

double shooting(double x0, double xn, double y0, double yn, double z0, double h) {

double z1 = z0, z2 = z0 + 1, y1, y2;

while (1) {

double y = y0, z = z1, x = x0;

for (; x < xn; x += h) rk4(&y, &z, x, h);

y1 = y;
 y = y0, z = z2, x = x0;

for (; x < xn; x += h) rk4(&y, &z, x, h);

y2 = y;

if (fabs(y1 - yn) < TOL) return z1;

z1 = z2 - (y2 - yn) * (z2 - z1) / (y2 - y1);

z2 += 1;

int main() {

double x0 = 0, xn = 1, y0 = 0, yn = 1, z0 = 1, h = 0.1;

double z = shooting(x0, xn, y0, yn, z0, h);

printf("Approximate slope at x0: %f\n", z);

return 0;

Python code:

import numpy as np

def f(x, y, z):

return z

def g(x, y, z):

return -y

def rk4(y, z, x, h):

k1, l1 = h * f(x, y, z), h * g(x, y, z)
k2, l2 = h * f(x + h/2, y + k1/2, z + l1/2), h * g(x + h/2, y + k1/2,
z + l1/2)
k3, l3 = h * f(x + h/2, y + k2/2, z + l2/2), h * g(x + h/2, y + k2/2,
z + l2/2)
k4, l4 = h * f(x + h, y + k3, z + l3), h * g(x + h, y + k3, z + l3)

y_new = y + (k1 + 2*k2 + 2*k3 + k4) / 6

z_new = z + (l1 + 2*l2 + 2*l3 + l4) / 6
 x_new = x + h

return y_new, z_new, x_new

def shooting(x0, xn, y0, yn, z0, h, tol=1e-6):

z1, z2 = z0, z0 + 1

while True:
y, z, x = y0, z1, x0
while x < xn:
y, z, x = rk4(y, z, x, h)
y1 = y

y, z, x = y0, z2, x0
while x < xn:
y, z, x = rk4(y, z, x, h)
y2 = y

if abs(y1 - yn) < tol:

return z1

z1, z2 = z2 - (y2 - yn) * (z2 - z1) / (y2 - y1), z2 + 1

x0, xn, y0, yn, z0, h = 0, 1, 0, 1, 1, 0.1

print(f"Approximate slope at x0: {shooting(x0, xn, y0, yn, z0, h)}")
R.K Method

C code:

#include <stdio.h>

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

void rungeKutta(float x0, float y0, float h, int n) {

for (int i = 0; i < n; i++) {

float k1 = h * f(x0, y0);

float k2 = h * f(x0 + h / 2, y0 + k1 / 2);

float k3 = h * f(x0 + h / 2, y0 + k2 / 2);

float k4 = h * f(x0 + h, y0 + k3);

y0 += (k1 + 2 * k2 + 2 * k3 + k4) / 6;

x0 += h;

printf("x = %.2f, y = %.4f\n", x0, y0);

int main() {

rungeKutta(0, 1, 0.1, 10);

return 0;

Python code:

def runge_kutta(x0, y0, h, n):

for i in range(n):
k1 = h * f(x0, y0)
k2 = h * f(x0 + h/2, y0 + k1/2)
k3 = h * f(x0 + h/2, y0 + k2/2)
k4 = h * f(x0 + h, y0 + k3)

y0 = y0 + (k1 + 2*k2 + 2*k3 + k4) / 6

x0 = x0 + h
 return y0

def f(x, y):

return -y

print(runge_kutta(0, 1, 0.1, 10))

Two Point Forward and Backward difference

C code:

#include <stdio.h>

double forward_diff(double f1, double f2, double h) { return (f2 - f1) / h; }

double backward_diff(double f1, double f2, double h) { return (f1 - f2) / h; }

int main() {

double x1, x2, f1, f2, h;

printf("Enter x1, f(x1), x2, f(x2): ");

scanf("%lf %lf %lf %lf", &x1, &f1, &x2, &f2);

h = x2 - x1;

printf("Forward Difference: %lf\n", forward_diff(f1, f2, h));

printf("Backward Difference: %lf\n", backward_diff(f2, f1, h));

return 0;

Python code:

def forward_diff(f1, f2, h): return (f2 - f1) / h

def backward_diff(f1, f2, h): return (f1 - f2) / h

x1, f1, x2, f2 = map(float, input("Enter x1, f(x1), x2, f(x2): ").split())

h = x2 - x1

print(f"Forward Difference: {forward_diff(f1, f2, h)}")

print(f"Backward Difference: {backward_diff(f2, f1, h)}")

Three Point Formula

C code:

#include <stdio.h>

double three_point_formula(double f_x1, double f_x2, double f_x3, double h) {

return (f_x1 - 4 * f_x2 + 3 * f_x3) / (2 * h);

int main() {

double f_x1, f_x2, f_x3, h;

printf("Enter f(x-2h), f(x-h), f(x) and step size h: ");

scanf("%lf %lf %lf %lf", &f_x1, &f_x2, &f_x3, &h);

printf("Derivative: %lf\n", three_point_formula(f_x1, f_x2, f_x3, h));

return 0;

Python code:

def three_point_formula(f_x1, f_x2, f_x3, h):

return (f_x1 - 4 * f_x2 + 3 * f_x3) / (2 * h)

f_x1, f_x2, f_x3, h = map(float, input("Enter f(x-2h), f(x-h), f(x), and

h: ").split())
print("Derivative:", three_point_formula(f_x1, f_x2, f_x3, h))
Derivatives using Newton’s Forward Difference

C code:

#include <stdio.h>

void newton_forward_derivative(float x[], float y[], int n, float h) {

for (int i = 0; i < n - 1; i++) y[i] = (y[i + 1] - y[i]) / h;

printf("First derivative at x0: %.5f\n", y[0]);

int main() {

int n;

printf("Enter number of data points: "); scanf("%d", &n);

float x[n], y[n], h;

printf("Enter x values: "); for (int i = 0; i < n; i++) scanf("%f", &x[i]);

printf("Enter y values: "); for (int i = 0; i < n; i++) scanf("%f", &y[i]);

h = x[1] - x[0];

newton_forward_derivative(x, y, n, h);

return 0;

Python code:

import numpy as np

def newton_forward_derivative(x, y):

return (y[1] - y[0]) / (x[1] - x[0])

x = np.array(list(map(float, input("Enter x values: ").split())))

y = np.array(list(map(float, input("Enter y values: ").split())))

print(f"First derivative at x0: {newton_forward_derivative(x, y):.5f}")
Derivatives using Newton’s backward difference

C code:

#include <stdio.h>

double newton_backward(double x[], double y[], int n, double h) {

for (int i = 1; i < n; i++)

for (int j = n - 1; j >= i; j--)

y[j] = (y[j] - y[j - 1]) / h;

return y[n - 1];

int main() {

int n;

printf("Enter number of data points: ");

scanf("%d", &n);

double x[n], y[n], h;

printf("Enter x values: ");

for (int i = 0; i < n; i++) scanf("%lf", &x[i]);

printf("Enter y values: ");

for (int i = 0; i < n; i++) scanf("%lf", &y[i]);

h = x[1] - x[0];

printf("Derivative at last point: %lf\n", newton_backward(x, y, n, h));

return 0;

Python code:

import numpy as np
def newton_backward(x, y):

n, h = len(y), x[1] - x[0]

for i in range(1, n):

y[n-1] = (y[n-1] - y[n-1-i]) / h

return y[-1]

n = int(input("Enter number of data points: "))

x = np.array(list(map(float, input("Enter x values: ").split())))

y = np.array(list(map(float, input("Enter y values: ").split())))

print("Derivative at last point:", newton_backward(x, y))

Romberg Intergration

C code:

#include <stdio.h>

#include <math.h>

#define f(x) (x * x)

double romberg(double a, double b, int n) {

double R[n][n], h = b - a;

R[0][0] = (f(a) + f(b)) * h / 2.0;

for (int i = 1; i < n; i++) {

h /= 2.0;

double sum = 0.0;

for (int k = 1; k <= (1 << (i - 1)); k++)

sum += f(a + (2 * k - 1) * h);

R[i][0] = 0.5 * R[i - 1][0] + h * sum;

for (int j = 1; j <= i; j++)

R[i][j] = (pow(4, j) * R[i][j - 1] - R[i - 1][j - 1]) / (pow(4, j) - 1);

return R[n - 1][n - 1];

int main() {

double a = 0, b = 1;

int n = 4;

printf("Romberg Integration: %.6f\n", romberg(a, b, n));

 return 0;

Python code:

import numpy as np

def f(x):

return x ** 2

def romberg(a, b, n):

R = np.zeros((n, n))

h=b-a

R[0, 0] = (f(a) + f(b)) * h / 2.0

for i in range(1, n):

h /= 2.0

sum_ = sum(f(a + (2 * k - 1) * h) for k in range(1, 2**(i-1) + 1))

R[i, 0] = 0.5 * R[i - 1, 0] + h * sum_

for j in range(1, i + 1):

R[i, j] = (4**j * R[i, j - 1] - R[i - 1, j - 1]) / (4**j - 1)

return R[n - 1, n - 1]

a, b, n = 0, 1, 4

print(f"Romberg Integration: {romberg(a, b, n):.6f}")

taylor’s seires:

c code:

#include <stdio.h>

double taylor_series(int x, int n) {

double res = 1, term = 1;

for (int i = 1; i < n; i++) {

term *= (double)x / i;

res += term;

return res;

int main() {

int x = 2, n = 10;

printf("e^%d ≈ %lf\n", x, taylor_series(x, n));

return 0;

Python code:

def taylor_series(x, n):

res, term = 1, 1

for i in range(1, n):

term *= x / i

res += term

return res

x, n = 2, 10

print(f"e^{x} ≈ {taylor_series(x, n)}")

Picard’s Method

C code:

#include <stdio.h>

#include <math.h>

#define f(x, y) (x + y)

void picard(double x0, double y0, double h, int n) {

double x = x0, y = y0;

printf("x = %.4f, y = %.4f\n", x, y);

for (int i = 0; i < n; i++) {

y = y0 + h * (x0 + y0);

x += h;

printf("x = %.4f, y = %.4f\n", x, y);

y0 = y;

int main() {

double x0 = 0, y0 = 1, h = 0.1;

int n = 5;

picard(x0, y0, h, n);

return 0;

Python code:

def picard(x0, y0, h, n):

x, y = x0, y0
 print(f"x = {x:.4f}, y = {y:.4f}")

for _ in range(n):

y = y0 + h * (x0 + y0)

x += h

print(f"x = {x:.4f}, y = {y:.4f}")

y0 = y

x0, y0, h, n = 0, 1, 0.1, 5

picard(x0, y0, h, n)
Laplace Method

C code:

#include <stdio.h>

#include <math.h>

#define N 10

#define EPS 1e-6

void solve_laplace(double grid[N][N]) {

int i, j;

double diff, maxDiff;

do {

maxDiff = 0.0;

for (i = 1; i < N - 1; i++) {

for (j = 1; j < N - 1; j++) {

double temp = grid[i][j];

grid[i][j] = 0.25 * (grid[i+1][j] + grid[i-1][j] + grid[i][j+1] + grid[i][j-1]);

diff = fabs(grid[i][j] - temp);

if (diff > maxDiff) maxDiff = diff;

} while (maxDiff > EPS);

void print_grid(double grid[N][N]) {

for (int i = 0; i < N; i++) {

for (int j = 0; j < N; j++) {

printf("%.2f ", grid[i][j]);

}
 printf("\n");

int main() {

double grid[N][N] = {0};

solve_laplace(grid);

print_grid(grid);

return 0;

Python code:

import numpy as np

N, EPS = 10, 1e-6

def solve_laplace(grid):

while True:

new_grid = grid.copy()

new_grid[1:-1, 1:-1] = 0.25 * (grid[:-2, 1:-1] + grid[2:, 1:-1] + grid[1:-1, :-2] + grid[1:-1, 2:])

if np.max(np.abs(new_grid - grid)) < EPS:

break

grid = new_grid

return grid

grid = np.zeros((N, N))

grid = solve_laplace(grid)

print(grid)
Poisson Method

C code:

#include <stdio.h>

#include <math.h>

#define N 50

#define MAX_ITER 10000

#define TOL 1e-6

void solve_poisson(double u[N][N], double f[N][N]) {

double u_new[N][N], error;

for (int iter = 0; iter < MAX_ITER; iter++) {

error = 0.0;

for (int i = 1; i < N - 1; i++) {

for (int j = 1; j < N - 1; j++) {

u_new[i][j] = 0.25 * (u[i+1][j] + u[i-1][j] + u[i][j+1] + u[i][j-1] - f[i][j]);

error += fabs(u_new[i][j] - u[i][j]);

if (error < TOL) break;

for (int i = 1; i < N - 1; i++)

for (int j = 1; j < N - 1; j++)

u[i][j] = u_new[i][j];

int main() {

double u[N][N] = {0}, f[N][N] = {0};

 f[N/2][N/2] = 1.0;

solve_poisson(u, f);

printf("Solution at a few points:\n");

for (int i = N/2 - 2; i <= N/2 + 2; i++) {

for (int j = N/2 - 2; j <= N/2 + 2; j++) {

printf("%6.3f ", u[i][j]);

printf("\n");

return 0;

Python code:

import numpy as np

N, MAX_ITER, TOL = 50, 10000, 1e-6

def solve_poisson():

u, f = np.zeros((N, N)), np.zeros((N, N))

f[N//2, N//2] = 1.0

for _ in range(MAX_ITER):

u_new = 0.25 * (np.roll(u, 1, axis=0) + np.roll(u, -1, axis=0) +

np.roll(u, 1, axis=1) + np.roll(u, -1, axis=1) - f)

error = np.abs(u_new - u).sum()

 if error < TOL:

break

u = u_new

return u

u = solve_poisson()

print("Solution at a few points:")

for i in range(N//2 - 2, N//2 + 3):

print(" ".join(f"{u[i, j]:6.3f}" for j in range(N//2 - 2, N//2 + 3)))

Cubic spline

C code:

#include <stdio.h>

#include <stdlib.h>

typedef struct {

double *a, *b, *c, *d, *x;

int n;

} Spline;

Spline* cubic_spline(double x[], double y[], int n) {

Spline *s = (Spline*)malloc(sizeof(Spline));

s->a = (double*)malloc(n * sizeof(double));

s->b = (double*)malloc(n * sizeof(double));

s->c = (double*)calloc(n, sizeof(double));

s->d = (double*)malloc(n * sizeof(double));

s->x = x;

s->n = n;

double *h = (double*)malloc((n - 1) * sizeof(double));

double *alpha = (double*)malloc(n * sizeof(double));

double *l = (double*)malloc(n * sizeof(double));

double *mu = (double*)malloc(n * sizeof(double));

double *z = (double*)malloc(n * sizeof(double));

for (int i = 0; i < n - 1; i++) {

h[i] = x[i + 1] - x[i];

s->a[i] = y[i];

}
 s->a[n - 1] = y[n - 1];

for (int i = 1; i < n - 1; i++)

alpha[i] = (3 / h[i]) * (y[i + 1] - y[i]) - (3 / h[i - 1]) * (y[i] - y[i - 1]);

l[0] = 1, mu[0] = 0, z[0] = 0;

for (int i = 1; i < n - 1; i++) {

l[i] = 2 * (x[i + 1] - x[i - 1]) - h[i - 1] * mu[i - 1];

mu[i] = h[i] / l[i];

z[i] = (alpha[i] - h[i - 1] * z[i - 1]) / l[i];

l[n - 1] = 1, z[n - 1] = 0, s->c[n - 1] = 0;

for (int j = n - 2; j >= 0; j--) {

s->c[j] = z[j] - mu[j] * s->c[j + 1];

s->b[j] = (y[j + 1] - y[j]) / h[j] - h[j] * (s->c[j + 1] + 2 * s->c[j]) / 3;

s->d[j] = (s->c[j + 1] - s->c[j]) / (3 * h[j]);

free(h), free(alpha), free(l), free(mu), free(z);

return s;

void free_spline(Spline *s) {

free(s->a), free(s->b), free(s->c), free(s->d), free(s);

int main() {

double x[] = {0, 1, 2, 3};

 double y[] = {1, 2, 0, 2};

int n = 4;

Spline *s = cubic_spline(x, y, n);

for (int i = 0; i < n - 1; i++)

printf("Spline segment %d: a=%.2f, b=%.2f, c=%.2f, d=%.2f\n",

i, s->a[i], s->b[i], s->c[i], s->d[i]);

free_spline(s);

return 0;

Python code:

import numpy as np

def cubic_spline(x, y):

n = len(x)

h = np.diff(x)

a, b, c, d = np.array(y, dtype=float), np.zeros(n), np.zeros(n), np.zeros(n)

alpha = np.zeros(n)

for i in range(1, n - 1):

alpha[i] = (3 / h[i]) * (y[i + 1] - y[i]) - (3 / h[i - 1]) * (y[i] - y[i - 1])

l, mu, z = np.ones(n), np.zeros(n), np.zeros(n)

for i in range(1, n - 1):

l[i] = 2 * (x[i + 1] - x[i - 1]) - h[i - 1] * mu[i - 1]

mu[i] = h[i] / l[i]

z[i] = (alpha[i] - h[i - 1] * z[i - 1]) / l[i]

 c[n - 1] = 0

for j in range(n - 2, -1, -1):

c[j] = z[j] - mu[j] * c[j + 1]

b[j] = (y[j + 1] - y[j]) / h[j] - h[j] * (c[j + 1] + 2 * c[j]) / 3

d[j] = (c[j + 1] - c[j]) / (3 * h[j])

return list(zip(a, b, c, d))

x = [0, 1, 2, 3]

y = [1, 2, 0, 2]

splines = cubic_spline(x, y)

for i, (a, b, c, d) in enumerate(splines[:-1]):

print(f"Spline segment {i}: a={a:.2f}, b={b:.2f}, c={c:.2f}, d={d:.2f}")