Question

Using MATLAB, Consider the polynomial f(x) = 3x^3 + 5x^2 − 58x − 40. Find the three roots of the polynomial, i.e, x where f(x) = 0, using: (i) Bisection method, and (ii) Newton’s method. Report the number of iterations taken by each algorithm using a tolerance of 10^−8 .

Answer 1

main_script.m

clc;clear;close all;
f = @(x) 3*x^3 + 5*x^2 - 58*x - 40; %function
df = @(x) 9*x^2 + 10*x - 58; %derivative
eps = 1e-8;
max_iter=100;
hold on;
fplot(f,[-10,10]); %graph to guess roots
fplot(@(x) 0, [-10,10]); %x axis
[root1, iter1] = newton(f,df,-6,eps,max_iter); %first root
[root2, iter2] = newton(f,df,0,eps,max_iter); %second root
[root3, iter3] = newton(f,df,6,eps,max_iter); %third root

fprintf("Newton Method\n");
fprintf("Root: %f, Iterations:%d\n",root1,iter1);
fprintf("Root: %f, Iterations:%d\n",root2,iter2);
fprintf("Root: %f, Iterations:%d\n",root3,iter3);

[root1, iter1] = bisection(f,-6,-3,eps,max_iter); %first root
[root2, iter2] = bisection(f,-3,1,eps,max_iter); %second root
[root3, iter3] = bisection(f,3,10,eps,max_iter); %third root

fprintf("Bisection Method\n");
fprintf("Root: %f, Iterations:%d\n",root1,iter1);
fprintf("Root: %f, Iterations:%d\n",root2,iter2);
fprintf("Root: %f, Iterations:%d\n",root3,iter3);

--------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------

newton.m

function [ root, iter ] = newton( f, df, x0, epsilon, max_iter )
%function, derivative,initial guess, epsilonerance, max iteration
x(1) = x0 - (f(x0)/df(x0));
error(1) = abs(x(1)-x0);
k = 2;
while (error(k-1) >= epsilon) && (k <= max_iter) %while we have a reasonable root
x(k) = x(k-1) - (f(x(k-1))/df(x(k-1))); %formula
error(k) = abs(x(k)-x(k-1)); %error
k = k+1;
end
root = x(length(x));
iter=k-2;
end

----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------

bisection.m

function [p,iter] = bisection(f,a,b, eps, max_iter)

prev=1000;
p = (a + b)/2;
err = 100;
iter=1;
while err >= eps && iter<max_iter
if f(a)*f(p)<0
b = p;
else
a = p;
end
prev=p;
p = (a + b)/2;
err = abs((prev-p));
%fprintf("%d Error:%.4f\n",iter,err);
iter = iter+1;
end

%fprintf("%d Error:%f\n",iter,err);
end