Effacer les filtres
Effacer les filtres

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

23 vues (au cours des 30 derniers jours)
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

Réponse acceptée

Torsten
Torsten le 18 Sep 2023
Modifié(e) : Torsten le 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
  1 commentaire
ABDALLA AL KHALEDI
ABDALLA AL KHALEDI le 19 Sep 2023
Modifié(e) : ABDALLA AL KHALEDI le 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;

Connectez-vous pour commenter.

Plus de réponses (1)

Sam Chak
Sam Chak le 18 Sep 2023
Modifié(e) : Sam Chak le 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
  2 commentaires
Dyuman Joshi
Dyuman Joshi le 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
ABDALLA AL KHALEDI le 19 Sep 2023
@Sam Chak thank you for your input

Connectez-vous pour commenter.

Produits


Version

R2019b

Community Treasure Hunt

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

Start Hunting!

Translated by