Numerical Method for Root Findings

Introduction to Polynomial Equations

Polynomial equations appear in various engineering and scientific problems, such as:
Modeling physical systems (e.g., motion equations in mechanics)
Electrical circuit analysis (e.g., polynomial equations for impedance)
Control systems (e.g., stability analysis using characteristic equations)
Optimization problems

A general polynomial equation in one variable can written as

where are coefficients.

Representing Polynomials
Polynomials are represented as vectors of coefficients ordered from the highest power to the constant term.

Example:

For , use:

p = [1, -6, 11, -6];

Finding Roots

Use roots() to find the roots of a polynomial.

 roots([1, -6, 11, -6])

ans = 3×1
3.0000
2.0000
1.0000

Differentiation

Compute the derivative using polyder()

dp = polyder([1, -6, 11, -6])

dp = 1×3
3 -12 11

[3, -12, 11] when converted to Polynomial equation based on its coefficient will become

that is equal to

Integration

Integrate a polynomial with polyint(p, k), where k is the constant of integration (default: 0).

Integrate: and take the constant of integration as 0.

ip = polyint([1, -6, 11, -6])

ip = 1×5
0.2500 -2.0000 5.5000 -6.0000 0

Polynomial Multiplication

Multiply polynomials using conv().

Multiply
 product = conv([1, 1], [1, -1])

product = 1×3
1 0 -1

Polynomial Division

Divide polynomials with deconv() to get quotient and remainder.

Example:

Divide by :

[quotient, remainder] = deconv([1, 0, -1], [1, -1])

quotient = 1×2
1 1
remainder = 1×3
0 0 0

Evaluating Polynomials

Evaluate at a point using polyval(p, x).

Evaluate at

x= [1, 3 ,5, 7, 10];

value = polyval(p,x)

value = 1×5
0 0 24 120 504

Polynomial giving 0 at 1 and 3 as they are the roots of this polynomial. How to find the roots, we will learn in the subsquent section.

Introduction to Root Finding

Root finding is a fundamental problem in numerical analysis, where the goal is to determine the values of variables that satisfy a given equation or system
of equations. These roots are the solutions to the equation for a single-variable function, or and for a system of equations
in multiple variables.

Solving Polynomial Equations Using the roots Function

In MATLAB, the roots function is used to find the roots of a single-variable polynomial. The function takes a vector of coefficients as input and returns the
roots of the polynomial.

Solve:

coefficients = [1, -7, 15.75, -11.25]; % Define the coefficients of the polynomial

roots_of_polynomial = roots(coefficients); % Find the roots using the roots function

% The roots of the polynomial are:

disp(roots_of_polynomial);

3.0000
2.5000
1.5000

This works well for polynomials of moderate degree. However, for high-degree polynomials or polynomials with complex coefficients, numerical stability
issues may arise.

NOTE: The roots function is ideal for solving single-variable polynomial equations. When roots fails to provide a solution for complex problems,
numerical methods such as the Bisection Method and Newton-Raphson Method are often employed. These methods are robust and can handle cases
where analytical or built-in functions like roots struggle.

Curve Fitting
Fit a polynomial to data using polyfit(x, y, n).

where x and y are the data points and n is order of polynomial.

 x = [1, 2, 3]; y = [5, 8, 11];
coeffs = polyfit(x, y, 1)

coeffs = 1×2
3.0000 2.0000

Returns [2, 3] meaning the above data can be represented as (y = 2x + 3)

Similarly, you can fit a very complex data as shown blow:

Let us generate some x and y data and then ask the polyfit to give us polynomial fit of n order. Once we visualize it we will know if the fitted equation is
good or not.
% 100 points between 0 and 2.5
x = linspace(0, 2.5, 100);

% Add random noise (0 to 0.1)

y = sin(x) + rand(1, length(x))/10;

% Fit 2nd-order polynomial to the data

p = polyfit(x, y, 2);

% Evaluate the polynomial using polyval

y_fit = polyval(p, x); % This is our equation which is being evaluated at every point of x.

figure;
plot(x, y, '-k', 'MarkerSize', 6, 'DisplayName', 'Noisy Data'); % Original data
hold on;
plot(x, y_fit, 'LineWidth', 2, 'DisplayName', '2nd-Order Fit'); % Fitted curve
hold off;
xlabel('x');
ylabel('y');
title('Noisy Sine Wave vs. Polynomial Fit');
legend('Location', 'best');
grid on;
Change the range, instead of going from 0 to 2.5 ... set these values to 0 to 5 and see that if the quadratic equation still able to fit those data.

Change the order (n value) to fit higher order polynomical equation.

Bisection Method
The bisection method is a root-finding technique that repeatedly bisects an interval and then selects a subinterval in which a root must lie for further
processing. Read more about it here: Bisection Method.

Following is the bare minimum code that you must know and it has been practiced in the class to solve differerent cases.
% Define the polynomial equation
func = @(x) x.^3 - x - 2;

% Initial interval [a, b]

a = 1;
b = 2;

% Tolerance and maximum iterations

tolerance = 1e-6; % or 0.000001
max_iter = 100;

% Initialize iteration counter

iteration = 0;

% Bisection algorithm
for k = 1:max_iter
iteration = iteration + 1;
c = (a + b) / 2; % Midpoint

% Check if midpoint is the exact root (unlikely in practice)

if func(c) == 0
break;
end

% Check if tolerance is met (exit condition)

if abs(c - a) < tolerance
break;
end

% Update the interval [a, b]

if func(a) * func(c) < 0
b = c; % Root lies in [a, c]
else
a = c; % Root lies in [c, b]
end
end

fprintf('Root found at x = %.6f after %d iterations.\n', c, iteration);

 Root found at x = 1.521380 after 20 iterations.

Its always a good way to visualize the function, if possible, that way you can find the interval for which the function has the opposite sign. To plot the above
function we can do this

f= @(x) x.^3 - x - 2;
xrange = -1:0.1:3;
plot(xrange , f(xrange),'--xk') % plotting these points [(x1,f(x1)), (x2, f(x2), ... , and so on)]
yline(0) % This is a horizontal line passing from y=0.

You can see here this line crossing y = 0 line in between 1 and 2. Therefore this could be an valid interval.
Using Bisection Method as Function
The following function implement the bisection method. Save it with the same name: bisection_method to find the roots of a function in one variable:
function root = bisection_method(func, a, b, tolerance, max_iteration)
% bisection_method: Finds the root of a function using the bisection method.
%
% Inputs:
% func: A function handle to the function whose root is to be found.
% a: The lower bound of the interval.
% b: The upper bound of the interval.
% tolerance: The tolerance for the root (stopping criterion).
% max_iteration: The maximum number of iterations allowed.
%
% Outputs:
% root: The approximate root of the function.

% Check if the function changes sign over the interval [a, b]

if func(a) * func(b) >= 0
error('The function must have opposite signs at the endpoints a and b.');
end

% Initialize variables
itereration = 0;
root = (a + b) / 2;

% While loop
while (b - a) / 2 > tolerance && itereration < max_iteration

% increment iteration and bisects the interval

itereration = itereration + 1;
root = (a + b) / 2;

% Check if the root is found

if func(root) == 0
break;
end

% Update the interval

if func(a) * func(root) < 0
b = root;
else
a = root;
 end
end

% Display the result

fprintf('Approximate root: %.10f\n', root);
fprintf('Number of iterations: %d\n', itereration);
end

Explanation
func: The function handle for which you want to find the root.

a, b: The endpoints of the interval [a,b] where the function is having opposite sign.

tolerance: The tolerance level, which determines how close the result should be to the actual root.

max_iteration: The maximum number of iterations to perform.

Notes
The function func must be continuous on the interval [a,b].

The interval [a,b] must contain a root, i.e., func(a) and func(b) must have opposite signs.

The method will stop when the interval size is smaller than the specified tolerance or when the maximum number of iterations is reached.

How to use this function

Suppose you want to find the root of the function in the interval with a tolerance of and a maximum of 100 iterations:

% Define the function

f = @(x) x.^3 - x - 2;

% Call the bisection method

root = bisection_method(f, 1, 2, 1e-6, 100);

Approximate root: 1.5213794708

Number of iterations: 19

Solve
Let us visualize it.

g= @(x) 4*x.^2-x.^3;
xrange = -1:0.1:5;
plot(xrange , g(xrange),'--or') % change f to g as now the function has been defined as g only.
yline(0)

You can see that it has two solutions one is from -1 to 1 and another solution lies in between 4 to 5.

[Actually this polynomial is having three solutions namely 0, 0 , 4 due to its highest power being 3]
 root = bisection_method(g, -10, 5, 1e-6, 100)

Approximate root: 3.9999991655

Number of iterations: 23
root = 4.0000

You may have seen the tables for different iterations at their respective wikipedia pages, we can also
implement here to print the different iterations in the tabular format.
Bisection Method
function results = bisection_method2(f, a, b, tol, max_iter)
% Bisection Method for finding roots of a function
%
% Inputs:
% f - Function handle @(x) f(x)
% a, b - Initial interval [a, b] where f(a)*f(b) < 0
% tol - Tolerance for stopping criterion
% max_iter - Maximum number of iterations
%
% Output:
% results - Table containing iteration, a_n, b_n, c_n, and f(c_n)

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

error('The function must have opposite signs at a and b.');
end

% Initialize table to store results

results = [];
iter = 0;

while (b - a) / 2 > tol && iter < max_iter

c = (a + b) / 2; % Midpoint
fc = f(c);

% Store the current iteration values

results = [results; iter + 1, a, b, c, fc];

if fc == 0 % Found exact root

break;
elseif f(a) * fc < 0
b = c;
else
a = c;
end
iter = iter + 1;
end

% Final midpoint and function value

c = (a + b) / 2;
 fc = f(c);
results = [results; iter + 1, a, b, c, fc];

% Convert to table for better readability

results = array2table(results, 'VariableNames', {'Iteration', 'a_n', 'b_n', 'c_n', 'f(c_n)'});
end

How to use it.

f = @(x) x^3 - x - 2; % Define function

a = 1; b = 2; % Initial interval
tol = 1e-6; % Tolerance
max_iter = 100; % Maximum iterations

results = bisection_method2(f, a, b, tol, max_iter);

disp(results);

Iteration a_n b_n c_n f(c_n)

_________ ______ ______ ______ ___________


1 1 2 1.5 -0.125
2 1.5 2 1.75 1.6094
3 1.5 1.75 1.625 0.66602
4 1.5 1.625 1.5625 0.2522
5 1.5 1.5625 1.5312 0.059113
6 1.5 1.5312 1.5156 -0.034054
7 1.5156 1.5312 1.5234 0.01225
8 1.5156 1.5234 1.5195 -0.010971
9 1.5195 1.5234 1.5215 0.00062218
10 1.5195 1.5215 1.5205 -0.0051789
11 1.5205 1.5215 1.521 -0.0022794
12 1.521 1.5215 1.5212 -0.00082891
13 1.5212 1.5215 1.5214 -0.00010343
14 1.5214 1.5215 1.5214 0.00025935
15 1.5214 1.5214 1.5214 7.7956e-05
16 1.5214 1.5214 1.5214 -1.2739e-05
17 1.5214 1.5214 1.5214 3.2608e-05 
18 1 5214 1 5214 1 5214 9 9343e 06
Newton Raphson Method
The Newton-Raphson method is an iterative root-finding algorithm that uses the derivative of a function to approximate its roots. Read it here Newton's
Method

This can be implement as [This you must know].

% Define the polynomial equation and its derivative

f = @(x) x.^3 - x - 2;
f_prime = @(x) 3*x.^2 - 1; % Derivative of f

% Initial guess
x0 = 1;

% Tolerance and maximum iterations

tolerance = 1e-6;
max_iter = 100;

% Initialize iteration counter

iteration = 0;

% Newton-Raphson algorithm
for k = 1:max_iter

iteration = iteration + 1;

x1 = x0 - func(x0)/f_prime(x0); % Update guess

if abs(x1 - x0) < tolerance

break;
end
x0 = x1; % Set new guess
end

% Print the result

fprintf('Root found at x = %.6f after %d iterations.\n', x1, iteration);

Root found at x = 1.521380 after 6 iterations.

Finding Square Root
Find the square root of 2:

Now the same problem has been translated into finding a roots of a equation
 % Function definition
f = @(x) x^2 - 2;
f_prime = @(x) 2*x;

% Initial guess
x0 = 1.5;

% Tolerance level
tolerance = 1e-6;

% Maximum number of iterations

max_iterations = 100;

% Newton-Raphson iteration
for i = 1:max_iterations
x1 = x0 - f(x0) / f_prime(x0);

% Check if the result is within the tolerance level

if abs(x1 - x0) < tolerance
break;
end

% Update x0 for the next iteration

x0 = x1;
end

% Display the result

fprintf('The approximate square root of 2 is: %.6f\n', x1);

The approximate square root of 2 is: 1.414214

Even you can minimize this code to find the square root of a. Use
 % Define the number we want the square root of
a = 2;

% Initial guess
x0 = 1;

% Iterate using the formula

for k = 1:100

x0 = (1/2) * (x0 + a/x0); % Apply the formula x_{n+1} = (1/2)(x_n + 2/x_n)

end

% Print the result

fprintf('Square root of 2 is approximately %.6f.\n', x0);

Square root of 2 is approximately 1.414214.

However, its always good to include the conditions. So that we can go fewer iteration and the program will break from the loop once we got satisfactory
result.

Newton Raphson Method can also be implemented as MATLAB function

function root = newton_raphson(func, dfunc, x0, tol, max_iter)
% newton_raphson: Finds the root of a function using the Newton-Raphson method.
%
% Inputs:
% func: A function handle to the function whose root is to be found.
% dfunc: A function handle to the derivative of the function.
% x0: The initial guess for the root.
% tol: The tolerance for the root (stopping criterion).
% max_iter: The maximum number of iterations allowed.
%
% Outputs:
% root: The approximate root of the function.

% Initialize variables
iter = 0;
root = x0;

% Newton-Raphson loop
while iter < max_iter
% Evaluate the function and its derivative at the current guess
f_val = func(root);
df_val = dfunc(root);

% Check if the derivative is zero to avoid division by zero

if df_val == 0
error('Derivative is zero. Newton-Raphson method cannot continue.');
end

% Update the root using the Newton-Raphson formula

root_new = root - f_val / df_val;

% Check for convergence

if abs(root_new - root) < tol
root = root_new;
break;
end

% Update the root for the next iteration

root = root_new;
 iter = iter + 1;
end

% Display the result

fprintf('Root found: %.6f\n', root);
fprintf('Number of iterations: %d\n', iter);
end

Suppose you want to find the root of the function with its derivative , starting from an initial guess , with a
tolerance of and a maximum of 100 iterations:

% Define the function and its derivative

f = @(x) x^3 - x - 2;
df = @(x) 3*x^2 - 1;

% Call the Newton-Raphson method

root = newton_raphson(f, df, 1, 1e-6, 100);

Root found: 1.521380

Number of iterations: 5

Explanation
1. func: The function handle for which you want to find the root.
2. dfunc: The function handle for the derivative of the function.
3. x0: The initial guess for the root.
4. tol: The tolerance level, which determines how close the result should be to the actual root.
5. max_iter: The maximum number of iterations to perform.

Notes
The Newton-Raphson method requires the derivative of the function. If the derivative is not available, you may need to use a numerical
approximation (e.g., finite differences).
The method may fail if the derivative is zero at any point or if the initial guess is not close enough to the root.
The method converges quickly if the initial guess is good, but it may diverge for poorly chosen initial guesses.

Newton Raphson Method

function results = newton_raphson_method2(f, df, x0, tol, max_iter)
% Newton-Raphson Method for finding roots of a function
%
% Inputs:
% f - Function handle @(x) f(x)
% df - Derivative of the function @(x) df(x)
% x0 - Initial guess
% tol - Tolerance for stopping criterion
% max_iter - Maximum number of iterations
%
% Output:
% results - Table containing iteration, x_n, and f(x_n)

% Initialize variables
iter = 0; % Counter for the number of iterations
x_n = x0; % Initialize the current approximation to the initial guess
results = []; % Initialize an empty array to store results

while iter < max_iter

fx = f(x_n); % Evaluate the function at the current approximation
dfx = df(x_n); % Evaluate the derivative of the function at the current approximation

% Check if the derivative is close to zero (to prevent division by zero)

if abs(dfx) < eps
error('Derivative near zero. The method may fail to converge.');
end

% Store the current iteration values in the results array

results = [results; iter + 1, x_n, fx];

% Check for convergence: if the function value is less than the tolerance, stop
if abs(fx) < tol
break;
end

% Update the approximation using the Newton-Raphson formula

x_n = x_n - fx / dfx;
iter = iter + 1; % Increment the iteration counter
end
 % Convert the results array into a table for better readability
results = array2table(results, 'VariableNames', {'Iteration', 'x_n', 'f(x_n)'});
end

Print The Table

f = @(x) x^3 - x - 2; % Define the function

df = @(x) 3*x^2 - 1; % Define its derivative
x0 = 1.5; % Initial guess
tol = 1e-6; % Tolerance
max_iter = 100; % Maximum iterations

results = newton_raphson_method2(f, df, x0, tol, max_iter);

disp(results);

Iteration x_n f(x_n)

_________ ______ __________

1 1.5 -0.125
2 1.5217 0.0021369
3 1.5214 5.8939e-07

Visualize the iteration

function results = newton_raphson_method3(f, df, x0, tol, max_iter)
% Newton-Raphson Method for finding roots of a function
%
% Inputs:
% f - Function handle @(x) f(x)
% df - Derivative of the function @(x) df(x)
% x0 - Initial guess
% tol - Tolerance for stopping criterion
% max_iter - Maximum number of iterations
%
% Output:
% results - Table containing iteration, x_n, and f(x_n)

% Initialize variables
iter = 0; % Counter for the number of iterations
x_n = x0; % Initialize the current approximation to the initial guess
iterations = []; % Initialize an array to store iteration numbers
x_values = []; % Initialize an array to store x_n values
fx_values = []; % Initialize an array to store f(x_n) values

while iter < max_iter

fx = f(x_n); % Evaluate the function at the current approximation
dfx = df(x_n); % Evaluate the derivative of the function at the current approximation

% Check if the derivative is close to zero (to prevent division by zero)

if abs(dfx) < eps
error('Derivative near zero. The method may fail to converge.');
end

% Store the current iteration values

iterations = [iterations; iter + 1];
x_values = [x_values; x_n];
fx_values = [fx_values; fx];

% Check for convergence: if the function value is less than the tolerance, stop
if abs(fx) < tol
break;
end
 % Update the approximation using the Newton-Raphson formula
x_n = x_n - fx / dfx;
iter = iter + 1; % Increment the iteration counter
end

% Create a table with iteration, x_n, and f(x_n)

results = table(iterations, x_values, fx_values, 'VariableNames', {'Iteration', 'x_n', 'f_x_n'});
end

Visualize how the x_n and f(x_n) are varying.

f = @(x) x^3 - x - 2; % Define the function
df = @(x) 3*x^2 - 1; % Define its derivative
x0 = 1; % Initial guess
tol = 1e-6; % Tolerance
max_iter = 100; % Maximum iterations
results = newton_raphson_method3(f, df, x0, tol, max_iter);

% Plot x_n vs Iteration

figure;
plot(results.Iteration, results.x_n, '-o');
xlabel('Iteration');
ylabel('x_n');
title('Convergence of x_n');
grid on;
% Plot f(x_n) vs Iteration
figure;
plot(results.Iteration, results.f_x_n, '-o');
xlabel('Iteration');
ylabel('f(x_n)');
title('Convergence of f(x_n)');
grid on;