How to use Bisection Program to locate Roots of given function
22 views (last 30 days)
Show older comments
The problem gives a differential equation and asks to find the roots using Bisection Method implimented into a MatLab function. I believe I have the correct program typed in MatLab for finding roots using bisection but I am struggling with how to input the given equation and find results. The problem is attached as a pdf and my code for the program is below: Any help is appriciated, Thank you!
function[root, fx, ez, iter] = bisect(func, x1, xu, es, maxit, varargin)
% bisect: root location zeros
%% input
%func = name of function
%x1, xu = lower and upper guesses
%es = desired relative error
%maxit = max iterations
%p1,p2 = additional parameters used by func
%%output
%root = real root
%fx = function value at root
%ea = approx. relative error
%iter = number of iterations
if nargin<3, error('at least 3 input arguments required'), end
test = func(x1,varargin{:})*func(xu,varargin{:});
if test>0, error('no sign change'), end
iter=0; xr=x1; ea=100;
while(1)
xrold=xr;
xr=(x1+xu)/2;
iter=iter+1;
if xr~=0, ea=abs((xr-xrold)/xr)*100;end
test=func(x1,varargin{:})*func(xr,varargin{:})
if test<0
xu=xr;
elseif test>0
x1=xr;
else
ea=0;
end
if ea<=es | iter>=maxit,break,end
end
root=xr; fx=func(xr, varargin{:})
7 Comments
James Tursa
on 5 Feb 2021
"... when the amount of liquid in the tank is at a maximun/minimum ..."
OK, now I see what you are trying to do.
Answers (2)
Jaime Varela
on 5 Feb 2021
Edited: Jaime Varela
on 5 Feb 2021
One thing I could see is that you are trying to check exactly when the function is zero
else
ea = 0
Numerically it's unlikely that an aritrary function will yield zero. You need to check when the function is approximately zero. A 11 line implementation of the bisection method is found here:
0 Comments
James Tursa
on 5 Feb 2021
Edited: James Tursa
on 5 Feb 2021
If the derivative depended only on the time t then your strategy for finding the maximum/minimum points would be a good one. Just finding the roots of the derivative function would get you what you want.
However, this approach probably isn't going to give you anything useful in your particular situation because your derivative depends on both the time t and the state y. You would be solving this equation:
0 = 3*(Q/A)*sin(t)^2 - (alpha/A)*(1+y)^1.5
Or, equivalently with Q = 600 and alpha = 120
sin(t)^2 = (1/15)*(1+y)^1.5
The thing is, there are an infinite number of solutions to the above equation. It is not clear from the drawing how far negative y can be, but for a large percentage of t values you can find a y value that satisfies the above equation. You wouldn't even need the bisection method for this because you can solve for the solution directly (i.e., given a t value you can solve directly for y immediately). So a large percentage of t values would have a corresponding y value that makes the derivative 0. But do these infinite number of solutions actually mean anything to you? Probably not. Whether your system actually gets into one of those infinite solution states will depend entirely on initial conditions. Of the infinite solutions available, only some of those states get achieved. The only thing the above equation tells you is that if at time t the tank level is a particular y value, then the tank level will be at a local maximum/minimum. What it doesn't tell you is if that particular state will ever be achieved, because actually getting into that state will depend on the initial conditions. It also doesn't tell you what the initial conditions have to be to get into a particular maximum/minimum state at a specific time t.
What you probably need to do instead is run simulations of this with specific initial conditions, e.g. with ode45( ), and then look for the maximum/minimum points in the numeric solutions. Were you given initial conditions for this problem? I.e., a starting y value?
4 Comments
James Tursa
on 8 Feb 2021
I would assume you simply pass a function handle that you are trying to drive to 0. E.g.,
func = @(t) sin(t).^2 - (1/15)*(1+1)^1.5
See Also
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!