Código en matlab para la bisección de Chapra 2023

 

Código en matlab para la bisección de Chapra 2023

Ejecutar:

>> fm=@(m,cd,t,v) sqrt(9.81*m/cd)*tanh(sqrt(9.81*cd/m)*t)-v

fm =

function_handle with value:

@(m,cd,t,v)sqrt(9.81*m/cd)*tanh(sqrt(9.81*cd/m)*t)-v

>> [mass fx ea iter]=mn_bisect(@(m) fm(m,0.25,4,36),40,200)

mass =

142.7377


fx =

4.6089e-07


ea =

5.3450e-05


iter =

21

Función, se debe guardar con el nombre mn_bisect.m:

function [root,fx,ea,iter]=mn_bisect(func,xl,xu,es,maxit,varargin)

% bisect: root location zeroes

% [root,fx,ea,iter]=bisect(func,xl,xu,es,maxit,p1,p2,...):

% uses bisection method to find the root of func

% input:

% func = name of function

% xl, xu = lower and upper guesses

% es = desired relative error (default = 0.0001%)

% maxit = maximum allowable iterations (default = 50)

% p1,p2,... = additional parameters used by func

% output:

% root = root estimate

% fx = function value at root estimate

% ea = approximate relative error (%)

% iter = number of iterations

if nargin<3,error('at least 3 input arguments required'),end

test = func(xl,varargin{:})*func(xu,varargin{:});

if test>0,error('no sign change'),end

if nargin<4 || isempty(es), es=0.0001;end

if nargin<5 || isempty(maxit), maxit=50;end

iter = 0; xr = xl; ea = 100;

while (1)

xrold = xr; xr = (xl + xu)/2;

iter = iter + 1;

if xr ~= 0,ea = abs((xr - xrold)/xr) * 100;end

test = func(xl,varargin{:})*func(xr,varargin{:});

if test < 0

xu = xr;

elseif test > 0

xl = xr;

else

ea = 0;

end

if ea <= es || iter >= maxit,break,end

end

root = xr; fx = func(xr, varargin{:});