# Application of any numerical root finding method (secant, bisection, etc.)

5 views (last 30 days)
ABDALLA AL KHALEDI on 18 Sep 2023
so I have this function that is a function of unknown parameter k, I need to find it using any numerical method. However, nothing seem to work for my allowable input values, there would seem to always be a combination of inputs that result in imaginary k values (physically impossible), no results (error for fzero), or the initial values do not give a change of sign, I have tried every thing including consulting chatGPT. Help is aapriciated.
Cs = 0;
d = 0.5;
g = 9.81;
%User definition of input parameters
%Possible values H = 0.05, 0.1, 0.2, 0.25
prompt = 'What is The Wave Height [m] ? ';
H = input(prompt);
%Possible values T = 0.6, 0.8, 1.5, 2, 2.5, 3, 5, 6
prompt = 'What is The Wave Period [s] ? ';
T = input(prompt);
%Function to be numerically solved for
f = @(k) sqrt(k/g)*Cs - ((2*pi)/(T*sqrt(g*k))) + (sqrt(tanh(k*d)))...
+(k*H/2)^2*((sqrt(tanh(k*d))*((2+7*(sech(2*k*d))^2)/(4*(1-(sech(2*k*d)))^2)))...
+(-0.5*sqrt(coth(k*d)))/(k*d))...
+(k*H/2)^4*((sqrt(tanh(k*d))*(4+32*(sech(2*k*d))...
-116*(sech(2*k*d))^2-400*(sech(2*k*d))^3-71*(sech(2*k*d))^4 ...
+146*(sech(2*k*d))^5)/(32*(1-(sech(2*k*d)))^5))+((sqrt(coth(k*d))*(2+4*(sech(2*k*d))...
+(sech(2*k*d))^2+2*(sech(2*k*d))^3)/(8*(1-(sech(2*k*d)))^3))/(k*d)));
%Initial values of k
k_zero = 4*(pi^2)/((T^2)*g);
k_one = 4*(pi^2)/((T^2)*g)*(coth((2*pi/T)*sqrt(d/g))^3/2)^2/3; %preferable one
I have tried several methods, like fzero as in writing:
% Use fzero to find the root
initial_guess = k_zero;
k = fzero(f, k_zero);
fprintf('Root k_hi = %.10f\n', k);
the secant method as in writing (after the above cod):
e = 10^-12;
n = 20;
for i=1:n
k_hi_two = (k_hi_zero*f(k_hi_one)-k_hi_one*f(k_hi_zero))/(f(k_hi_one)-f(k_hi_zero));
sprintf('k_hi%d = %.10f\n',i,k_hi_two);
if abs((k_hi_two-k_hi_one)/k_hi_two) < e
break
end
k_hi_zero = k_hi_one;
k_hi_one = k_hi_two;
end
k_hi = k_hi_two;
and the bisection method as in writing an m file function
function [k,e] = bis(f,a,b,tol)
% function [k e] = bisect(f,a,b,tol)
% Performs bisection until the error is less than tol
% Inputs:
% f: a function
% a, b: left and right edges of the interval
% tol: tolerance for stopping (e.g., 1e-6)
% Outputs:
% k: the estimated solution of f(k) = 0
% e: an upper bound on the error
% evaluate at the ends and make sure there is a sign change
c = f(a);
d = f(b);
if c*d > 0
error('Function has the same sign at both endpoints.');
end
while (b - a) / 2 > tol
% find the middle and evaluate there
k = (a + b) / 2;
y = f(k);
if y == 0 % solved the equation exactly
a = k;
b = k;
break; % exits the while loop
end
% decide which half to keep, so that the signs at the ends differ
if c * y < 0
b = k;
else
a = k;
end
end
% set the best estimate for k and the error bound
k = (a + b) / 2;
e = (b - a) / 2;
end

Torsten on 18 Sep 2023
Edited: Torsten on 18 Sep 2023
Plot first, then solve by using the approximate root as initial guess:
Cs = 0;
d = 0.5;
g = 9.81;
%User definition of input parameters
%Possible values H = 0.05, 0.1, 0.2, 0.25
H = 0.2;
%Possible values T = 0.6, 0.8, 1.5, 2, 2.5, 3, 5, 6
T = 2;
%Function to be numerically solved for
f = @(k) sqrt(k/g)*Cs - ((2*pi)/(T*sqrt(g*k))) + (sqrt(tanh(k*d)))...
+(k*H/2)^2*((sqrt(tanh(k*d))*((2+7*(sech(2*k*d))^2)/(4*(1-(sech(2*k*d)))^2)))...
+(-0.5*sqrt(coth(k*d)))/(k*d))...
+(k*H/2)^4*((sqrt(tanh(k*d))*(4+32*(sech(2*k*d))...
-116*(sech(2*k*d))^2-400*(sech(2*k*d))^3-71*(sech(2*k*d))^4 ...
+146*(sech(2*k*d))^5)/(32*(1-(sech(2*k*d)))^5))+((sqrt(coth(k*d))*(2+4*(sech(2*k*d))...
+(sech(2*k*d))^2+2*(sech(2*k*d))^3)/(8*(1-(sech(2*k*d)))^3))/(k*d)));
k = 0.4:0.01:5;
plot(k,arrayfun(@(k)f(k),k)) k_zero = 1.0;
% Use fzero to find the root
initial_guess = 1;
k = fzero(f, initial_guess);
fprintf('Root k_hi = %.10f\n', k);
Root k_hi = 1.5026507884
ABDALLA AL KHALEDI on 19 Sep 2023
Edited: ABDALLA AL KHALEDI on 19 Sep 2023
thank you @Torsten your solution served as a graphical reference to my issue, the idea was that my initial values do not necessarily bracket the root rendering the bisection method useless, and then I learened that when given two initial values, the fzero function automatically assumes that they bracket the root. What I did was using a code for the secant method found from this video, developed by "ATTIQ IQBAL" the resulting code was as shown below, when the initial guess is carefully selected from your graphical method, the results were a match for several input combinations.
Cs = 0;
d = 0.5;
g = 9.81;
%User definition of input parameters
%Possible values H = 0.05, 0.1, 0.2, 0.25
%Possible values T = 0.6, 0.8, 1.5, 2, 2.5, 3, 5, 6
prompt = 'What is The Wave Height [m] ? ';
H = input(prompt);
prompt = 'What is The Wave Period [s] ? ';
T = input(prompt);
%hi = higher order
f = @(k_hi) sqrt(k_hi/g)*Cs - ((2*pi)/(T*sqrt(g*k_hi))) + (sqrt(tanh(k_hi*d)))...
+(k_hi*H/2)^2*((sqrt(tanh(k_hi*d))*((2+7*(sech(2*k_hi*d))^2)/(4*(1-(sech(2*k_hi*d)))^2)))...
+(-0.5*sqrt(coth(k_hi*d)))/(k_hi*d))...
+(k_hi*H/2)^4*((sqrt(tanh(k_hi*d))*(4+32*(sech(2*k_hi*d))...
-116*(sech(2*k_hi*d))^2-400*(sech(2*k_hi*d))^3-71*(sech(2*k_hi*d))^4 ...
+146*(sech(2*k_hi*d))^5)/(32*(1-(sech(2*k_hi*d)))^5))+((sqrt(coth(k_hi*d))*(2+4*(sech(2*k_hi*d))...
+(sech(2*k_hi*d))^2+2*(sech(2*k_hi*d))^3)/(8*(1-(sech(2*k_hi*d)))^3))/(k_hi*d)));
%Linear Approximation (Exact for Deep and Shallow Waters (at most 1.7% error in between, Fenton)
%Current and nonlinearities are neglected
k_hi_zero = 4*(pi^2)/((T^2)*g)*(coth((2*pi/T)*sqrt(d/g))^3/2)^2/3;
%Linear Approximation Deep Water
k_hi_one = 4*(pi^2)/((T^2)*g);
e = 10^-14;
n = 20;
for i=1:n
k_hi_two = (k_hi_zero*f(k_hi_one)-k_hi_one*f(k_hi_zero))/(f(k_hi_one)-f(k_hi_zero));
sprintf('k_hi%d = %.10f\n',i,k_hi_two)
if abs((k_hi_two-k_hi_one)/k_hi_two) < e
break
end
k_hi_zero = k_hi_one;
k_hi_one = k_hi_two;
end
k_hi = k_hi_two;

Sam Chak on 18 Sep 2023
Edited: Sam Chak on 18 Sep 2023
Take pride in the fact that your Bisection method yields the same result as fzero().
Cs = 0;
d = 0.5;
g = 9.81;
% User definition of input parameters
% Possible values H = 0.05, 0.1, 0.2, 0.25
H = 0.2;
% Possible values T = 0.6, 0.8, 1.5, 2, 2.5, 3, 5, 6
T = 2;
% Function to be numerically solved for
f = @(k) sqrt(k/g)*Cs - ((2*pi)/(T*sqrt(g*k))) + (sqrt(tanh(k*d)))...
+(k*H/2)^2*((sqrt(tanh(k*d))*((2+7*(sech(2*k*d))^2)/(4*(1-(sech(2*k*d)))^2)))...
+(-0.5*sqrt(coth(k*d)))/(k*d))...
+(k*H/2)^4*((sqrt(tanh(k*d))*(4+32*(sech(2*k*d))...
-116*(sech(2*k*d))^2-400*(sech(2*k*d))^3-71*(sech(2*k*d))^4 ...
+146*(sech(2*k*d))^5)/(32*(1-(sech(2*k*d)))^5))+((sqrt(coth(k*d))*(2+4*(sech(2*k*d))...
+(sech(2*k*d))^2+2*(sech(2*k*d))^3)/(8*(1-(sech(2*k*d)))^3))/(k*d)));
k = 0.4:0.01:5;
plot(k, arrayfun(@(k) f(k), k)), grid on
axis([0 3 -1 1])
title('Bisection method')
xline(1, '-.', 'a = 1')
xline(2, '-.', 'b = 2')
xlabel('k') % Use Bisection to find the root
a = 1;
b = 2;
tol = 1e-10;
[k, e] = bis(f, a, b, tol)
k = 1.5027
e = 5.8208e-11
fprintf('Root k_hi = %.10f\n', k);
Root k_hi = 1.5026507885
function [k,e] = bis(f,a,b,tol)
% function [k e] = bisect(f,a,b,tol)
% Performs bisection until the error is less than tol
% Inputs:
% f: a function
% a, b: left and right edges of the interval
% tol: tolerance for stopping (e.g., 1e-6)
% Outputs:
% k: the estimated solution of f(k) = 0
% e: an upper bound on the error
% evaluate at the ends and make sure there is a sign change
c = f(a);
d = f(b);
if c*d > 0
error('Function has the same sign at both endpoints.');
end
while (b - a) / 2 > tol
% find the middle and evaluate there
k = (a + b) / 2;
y = f(k);
if y == 0 % solved the equation exactly
a = k;
b = k;
break; % exits the while loop
end
% decide which half to keep, so that the signs at the ends differ
if c * y < 0
b = k;
else
a = k;
end
end
% set the best estimate for k and the error bound
k = (a + b) / 2;
e = (b - a) / 2;
end
Dyuman Joshi on 18 Sep 2023
@Sam Chak, fzero uses a combination of bisection, secant and some other methods, so it would be suprising if OP's bisection method didn't produce the same (or a similar) result as fzero's result.
ABDALLA AL KHALEDI on 19 Sep 2023
@Sam Chak thank you for your input

### Categories

Find more on Solver Outputs and Iterative Display in Help Center and File Exchange

R2019b

### Community Treasure Hunt

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

Start Hunting!