0% found this document useful (0 votes)
2 views24 pages

Numerical Methods for Root Finding

Computational code

Uploaded by

josephkurianup
Copyright
© All Rights Reserved
We take content rights seriously. If you suspect this is your content, claim it here.
Available Formats
Download as PDF, TXT or read online on Scribd
0% found this document useful (0 votes)
2 views24 pages

Numerical Methods for Root Finding

Computational code

Uploaded by

josephkurianup
Copyright
© All Rights Reserved
We take content rights seriously. If you suspect this is your content, claim it here.
Available Formats
Download as PDF, TXT or read online on Scribd

 Bisection method

#include <iostream>

#include <cmath>

using namespace std;

double f(double x) {

return x*x*x - x - 2;

double bisection(double a, double b, double tolerance) {

if (f(a) * f(b) >= 0) {

cout << "Invalid interval: f(a) and f(b) must have opposite signs!" << endl;

return -1;

double mid;

while ((b - a) >= tolerance) {

mid = (a + b) / 2.0;

if (f(mid) == 0.0)

break;

if (f(a) * f(mid) < 0)

b = mid;

else

a = mid;

}
return mid;

int main() {

double a = 1, b = 2;

double tolerance = 0.0001;

double root = bisection(a, b, tolerance);

if (root != -1) {

cout << "Approximate root = " << root << endl;

cout << "f(root) = " << f(root) << endl;

return 0;

 REGULA FALSI METHOD

#include <iostream>

#include <cmath>

using namespace std;

double f(double x) {

return x*x*x - x - 2;

double regulaFalsi(double a, double b, double tolerance, int maxIter) {


double fa = f(a);

double fb = f(b);

if (fa * fb >= 0) {

cout << "Invalid interval: f(a) and f(b) must have opposite signs!" << endl;

return -1;

double c = a, fc;

for (int iter = 1; iter <= maxIter; ++iter) {

if (fb - fa == 0) {

cout << "Division by zero in formula. Stopping." << endl;

return -1;

c = (a * fb - b * fa) / (fb - fa);

fc = f(c);

cout << "Iter " << iter << " | a=" << a << " b=" << b

<< " c=" << c << " f(c)=" << fc << endl;

if (fabs(fc) < tolerance || fabs(b - a) < tolerance) {

cout << "Converged after " << iter << " iterations." << endl;

return c;

if (fa * fc < 0) {

b = c;

fb = fc;
} else {

a = c;

fa = fc;

cout << "Max iterations reached. Returning best c." << endl;

return c;

int main() {

double a = 1.0, b = 2.0;

double tolerance = 1e-6;

int maxIter = 100;

double root = regulaFalsi(a, b, tolerance, maxIter);

if (root != -1) {

cout << "\nApproximate root = " << root << endl;

cout << "f(root) = " << f(root) << endl;

return 0;

NEWTON RAPHSON METHOD

#include <iostream>

#include <cmath>
using namespace std;

double f(double x) {

return x*x*x - x - 2;

double df(double x) {

return 3*x*x - 1;

double newtonRaphson(double x0, double tolerance) {

double x1;

while (true) {

double fx = f(x0);

double dfx = df(x0);

if (dfx == 0) {

cout << "Derivative is zero. Cannot continue." << endl;

return x0;

x1 = x0 - fx / dfx;

if (fabs(x1 - x0) < tolerance)

break;

x0 = x1;
}

return x1;

int main() {

double x0, tolerance;

cout << "Enter initial guess: ";

cin >> x0;

cout << "Enter tolerance (e.g. 0.0001): ";

cin >> tolerance;

double root = newtonRaphson(x0, tolerance);

cout << "Approximate root = " << root << endl;

cout << "f(root) = " << f(root) << endl;

return 0;

TRAPEZOIDAL METHOD

#include <iostream>

#include <cmath>

#include <iomanip>

using namespace std;

double f(double x) {
return x * x; // Example: f(x) = x^2

int main() {

double a, b;

int n;

cout << "Enter lower limit a: ";

cin >> a;

cout << "Enter upper limit b: ";

cin >> b;

cout << "Enter number of subintervals n: ";

cin >> n;

double h = (b - a) / n;

double sum = f(a) + f(b);

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

double x = a + i * h;

sum += 2 * f(x);

double result = (h / 2) * sum;

cout << fixed << setprecision(6);

cout << "\nApproximate integral = " << result << endl;

return 0;
}

Simpsons rule [1/3]

#include <iostream>

#include <cmath>

#include <iomanip>

using namespace std;

double f(double x) {

return x * x; // Example: f(x) = x^2

int main() {

double a, b;

int n;

cout << "Enter lower limit a: ";

cin >> a;

cout << "Enter upper limit b: ";

cin >> b;

cout << "Enter number of subintervals n (even number): ";

cin >> n;

if (n % 2 != 0) {

cout << "Number of subintervals must be even for Simpson's 1/3 rule." << endl;

return 1;
}

double h = (b - a) / n;

double sum = f(a) + f(b);

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

double x = a + i * h;

if (i % 2 == 0) {

sum += 2 * f(x);

} else {

sum += 4 * f(x);

double integral = (h / 3) * sum;

cout << fixed << setprecision(10);

cout << "\nApproximate integral = " << integral << endl;

return 0;

Gauss quadrature method

#include <iostream>
#include <cmath>
using namespace std;

double f(double x) {
return x * x;
}

int main() {
double a, b;
cout << "Enter lower limit a: ";
cin >> a;
cout << "Enter upper limit b: ";
cin >> b;

double w1 = 1.0, w2 = 1.0;


double x1 = -1.0 / sqrt(3);
double x2 = 1.0 / sqrt(3);

double xm = 0.5 * (b + a);


double xr = 0.5 * (b - a);

double integral = xr * (w1 * f(xm + xr * x1) + w2 * f(xm + xr * x2));

cout << "Approximate integral using 2-point Gauss Quadrature = " << integral
<< endl;

return 0;
}

Gauss seidel method

#include <iostream>

#include <cmath>

using namespace std;

int main() {

int n;

cout << "Enter number of unknowns: ";


cin >> n;

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

cout << "Enter the coe icients of the matrix A row-wise:" << endl;

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

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

cin >> A[i][j];

cout << "Enter the constants vector b:" << endl;

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

cin >> b[i];

x[i] = 0;

double tol;

cout << "Enter tolerance (e.g. 0.0001): ";

cin >> tol;

int maxIter = 100;

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

double maxDi = 0;

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

double sum = 0;

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

if (j != i)

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


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

maxDi = max(maxDi , fabs(x_new - x[i]));

x[i] = x_new;

if (maxDi < tol)

break;

cout << "Solution:" << endl;

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

cout << "x" << i + 1 << " = " << x[i] << endl;

return 0;

Guass elimination

#include <iostream>

#include <cmath>

using namespace std;

int main() {

int n;

cout << "Enter number of unknowns: ";

cin >> n;
double A[10][10], b[10], x[10];

cout << "Enter the coe icients of matrix A row-wise:" << endl;

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

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

cin >> A[i][j];

cout << "Enter constants vector b:" << endl;

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

cin >> b[i];

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

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

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

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

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

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

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

x[i] = b[i];

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

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

x[i] /= A[i][i];

}
cout << "Solution:" << endl;

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

cout << "x" << i + 1 << " = " << x[i] << endl;

return 0;

Runga kutta method

#include <iostream>

#include <iomanip>

using namespace std;

double f(double x, double y) {

return x + y; // Example: dy/dx = x + y

int main() {

double x0, y0, xn, h;

cout << "Enter initial x (x0): ";

cin >> x0;

cout << "Enter initial y (y0): ";

cin >> y0;

cout << "Enter final x (xn): ";

cin >> xn;


cout << "Enter step size h: ";

cin >> h;

double n = (xn - x0) / h;

double x = x0, y = y0;

cout << fixed << setprecision(6);

cout << "x\t\ty" << endl;

cout << x << "\t" << y << endl;

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

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

double k2 = h * f(x + h/2.0, y + k1/2.0);

double k3 = h * f(x + h/2.0, y + k2/2.0);

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

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

x = x + h;

cout << x << "\t" << y << endl;

return 0;

Tylour series method


#include <iostream>

#include <iomanip>

using namespace std;

double f(double x, double y) {

return x + y;

double f_prime(double x, double y) {

return 1 + f(x, y);

int main() {

double x0, y0, xn, h;

cout << "Enter initial x (x0): ";

cin >> x0;

cout << "Enter initial y (y0): ";

cin >> y0;

cout << "Enter final x (xn): ";

cin >> xn;

cout << "Enter step size h: ";

cin >> h;

double n = (xn - x0) / h;

double x = x0, y = y0;

cout << fixed << setprecision(6);


cout << "x\t\ty" << endl;

cout << x << "\t" << y << endl;

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

double y_new = y + h * f(x, y) + (h*h/2.0) * f_prime(x, y);

x = x + h;

y = y_new;

cout << x << "\t" << y << endl;

return 0;

Matrix multiplication

#include <iostream>

using namespace std;

int main() {

int r1, c1, r2, c2;

cout << "Enter rows and columns of first matrix: ";

cin >> r1 >> c1;

cout << "Enter rows and columns of second matrix: ";

cin >> r2 >> c2;


if (c1 != r2) {

cout << "Matrix multiplication not possible (columns of first must equal rows of
second)." << endl;

return 0;

int A[10][10], B[10][10], C[10][10] = {0};

cout << "Enter elements of first matrix:" << endl;

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

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

cin >> A[i][j];

cout << "Enter elements of second matrix:" << endl;

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

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

cin >> B[i][j];

// Matrix multiplication

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

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

C[i][j] = 0;

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

C[i][j] += A[i][k] * B[k][j];

}
cout << "Resultant matrix:" << endl;

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

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

cout << C[i][j] << " ";

cout << endl;

return 0;

Matrix addition

#include <iostream>

using namespace std;

int main() {

int r, c;

cout << "Enter number of rows and columns: ";

cin >> r >> c;

int A[10][10], B[10][10], C[10][10];

cout << "Enter elements of first matrix:" << endl;

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

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

cin >> A[i][j];


cout << "Enter elements of second matrix:" << endl;

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

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

cin >> B[i][j];

// Matrix addition

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

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

C[i][j] = A[i][j] + B[i][j];

cout << "Resultant matrix after addition:" << endl;

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

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

cout << C[i][j] << " ";

cout << endl;

return 0;

Transpose

#include <iostream>

using namespace std;


int main() {

int r, c;

cout << "Enter number of rows and columns: ";

cin >> r >> c;

int A[10][10], T[10][10];

cout << "Enter elements of the matrix:" << endl;

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

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

cin >> A[i][j];

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

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

T[j][i] = A[i][j];

cout << "Transpose of the matrix:" << endl;

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

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

cout << T[i][j] << " ";

cout << endl;

return 0;

Substraction
#include <iostream>

using namespace std;

int main() {

int r, c;

cout << "Enter number of rows and columns: ";

cin >> r >> c;

int A[10][10], B[10][10], C[10][10];

cout << "Enter elements of first matrix:" << endl;

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

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

cin >> A[i][j];

cout << "Enter elements of second matrix:" << endl;

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

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

cin >> B[i][j];

// Matrix subtraction

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

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

C[i][j] = A[i][j] - B[i][j];

cout << "Resultant matrix after subtraction:" << endl;

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

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


cout << C[i][j] << " ";

cout << endl;

return 0;

Vector (addition, dot product cross product)

#include <iostream>

using namespace std;

int main() {

int n;

cout << "Enter dimension of vectors: ";

cin >> n;

double A[10], B[10], sum[10];

cout << "Enter elements of first vector:" << endl;

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

cin >> A[i];

cout << "Enter elements of second vector:" << endl;

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

cin >> B[i];

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

sum[i] = A[i] + B[i];


cout << "Vector addition result: ";

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

cout << sum[i] << " ";

cout << endl;

double dot = 0;

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

dot += A[i] * B[i];

cout << "Dot product: " << dot << endl;

if (n == 3) {

double cross[3];

cross[0] = A[1]*B[2] - A[2]*B[1];

cross[1] = A[2]*B[0] - A[0]*B[2];

cross[2] = A[0]*B[1] - A[1]*B[0];

cout << "Cross product: ";

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

cout << cross[i] << " ";

cout << endl;

} else {

cout << "Cross product is only defined for 3D vectors." << endl;

return 0;

You might also like