How to add the count of iterations and relative error to the code?

9 views (last 30 days)
Lingbai Ren
Lingbai Ren on 21 Oct 2021
I have this Newton-Rhapson and bisection method functions, but they cannot track how mant literations operated and what the relative error is. I am wondering where I can add lines of code to let them able to achieve those?
function root = newton(fun,x0,fprime,tol,maxiter)
% newton function computes the root of the function
% fun using Newton-Raphson method
% Check whether the initial guess is a root
% Error occurs if root is not determined within the
% desired tolerance in the max number iterations
root = x0;
iter = 0;
if ~isa(fprime, 'function_handle')
fprime = @(x) central_difference(fun,x,1e-4);
end
while iter < maxiter
if fprime(root) == x0
disp('The derivative is equal to 0, the current root is %.2f\n',root)
break
end
dx = -fun(root)/fprime(root);
root = root + dx;
if abs(dx) <= tol
break
end
iter = iter + 1;
end
if iter == maxiter && abs(dx) > tol
error('root not found in given iterations')
end
end
function [c] = bisection(fun,a,b,tol,maxlits,plot_output)
% bisection function computes the root, c, of a function, fun
% initial intercal
c_old = 0;
iters = 0;
e_a_vec = [];
% While loop
while iters < maxlits
iters = iters + 1;
c = (a + b)/2;
if plot_output == 1
x = linspace(a,b);
y = [];
for i = 1:length(x)
y(i) = fun(x(i));
end
nexttile
hold on
plot(a,fun(a),'-o')
plot(b,fun(b),'-o')
plot(c,fun(c),'-o')
plot(x,y,'k-')
yline(0,'--')
end
if iters >= 1
e_a = (c - c_old)/c;
e_a_vec(iters) = e_a;
if abs(b-a)/2 < tol | abs(fun(c)) <= tol | abs(e_a) <= tol
break
end
else abs(fun(c)) > tol
error('The root is not within the desired tolerance')
end
if fun(a) * fun(b) < 0
if fun(a) * fun(c) > 0
a = c;
else
b = c;
end
c_old = c_old + c;
else
error('The limites do not bracket a root')
end
end

Answers (1)

Sulaymon Eshkabilov
Sulaymon Eshkabilov on 21 Oct 2021
Here are how you can obtain the final iteration and rel./abs error values when the roots are found.
function root = newton(fun,x0,fprime,tol,maxiter)
% newton function computes the root of the function
% fun using Newton-Raphson method
% Check whether the initial guess is a root
% Error occurs if root is not determined within the
% desired tolerance in the max number iterations
root = x0;
iter = 0;
if ~isa(fprime, 'function_handle')
fprime = @(x) central_difference(fun,x,1e-4);
end
while iter < maxiter
if fprime(root) == x0
disp('The derivative is equal to 0, the current root is %.2f\n',root)
break
end
dx = -fun(root)/fprime(root);
root = root + dx;
if abs(dx) <= tol
fprintf('Simulation is halted after %d iterations \n', iter) % Number of iterations
fprintf('Simulation is halted after %f of Rel. Error (Derivative Value) \n', dx) % Rel. Error
break
end
iter = iter + 1;
end
if iter == maxiter && abs(dx) > tol
error('root not found in given iterations')
end
end
function [c] = bisection(fun,a,b,tol,maxlits,plot_output)
% bisection function computes the root, c, of a function, fun
% initial intercal
c_old = 0;
iters = 0;
e_a_vec = [];
% While loop
while iters < maxlits
iters = iters + 1;
c = (a + b)/2;
if plot_output == 1
x = linspace(a,b);
y = [];
for i = 1:length(x)
y(i) = fun(x(i));
end
nexttile
hold on
plot(a,fun(a),'-o')
plot(b,fun(b),'-o')
plot(c,fun(c),'-o')
plot(x,y,'k-')
yline(0,'--')
end
if iters >= 1
e_a = (c - c_old)/c;
e_a_vec(iters) = e_a;
if abs(b-a)/2 < tol | abs(fun(c)) <= tol | abs(e_a) <= tol
fprintf('Simulation is halted after %d iterations \n', iters) % Number of iterations
fprintf('Simulation is halted after %f of Rel Diff \n', abs(b-a)/2) % Rel. Difference
fprintf('Simulation is halted after %f of Abs. Fcn. Value \n', abs(fun(c))) % Abs. Value
fprintf('Simulation is halted after %f of Abs. Error \n', abs(e_a)) % Abs. Error
break
end
else abs(fun(c)) > tol
error('The root is not within the desired tolerance')
end
if fun(a) * fun(b) < 0
if fun(a) * fun(c) > 0
a = c;
else
b = c;
end
c_old = c_old + c;
else
error('The limites do not bracket a root')
end
end

Categories

Products


Release

R2021a

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!