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

조회 수: 7 (최근 30일)
Lingbai Ren
Lingbai Ren 2021년 10월 21일
답변: Sulaymon Eshkabilov 2021년 10월 21일
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

답변 (1개)

Sulaymon Eshkabilov
Sulaymon Eshkabilov 2021년 10월 21일
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

카테고리

Help CenterFile Exchange에서 Symbolic Math Toolbox에 대해 자세히 알아보기

제품


릴리스

R2021a

Community Treasure Hunt

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

Start Hunting!

Translated by