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

18 次查看（过去 30 天）
ABDALLA AL KHALEDI
ABDALLA AL KHALEDI 2023-9-18
评论： ABDALLA AL KHALEDI 2023-9-19
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
Torsten 2023-9-18
编辑：Torsten 2023-9-18
Ran in:
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
  1 个评论
ABDALLA AL KHALEDI
ABDALLA AL KHALEDI 2023-9-19
编辑：ABDALLA AL KHALEDI 2023-9-19
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;

请先登录，再进行评论。

更多回答（1 个）

Sam Chak
Sam Chak 2023-9-18
编辑：Sam Chak 2023-9-18
Ran in:
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
  2 个评论
Dyuman Joshi
Dyuman Joshi 2023-9-18
@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.

请先登录，再进行评论。

类别

在 Help Center 和 File Exchange 中查找有关 Calculus 的更多信息

标签

产品


版本

R2019b

Community Treasure Hunt

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

Start Hunting!

Translated by