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

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

Answer 1

[Torsten]

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]

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;

Answer 2

[Sam Chak]

Hi @ABDALLA AL KHALEDI

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]

@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]

@Sam Chak thank you for your input