Application of any numerical root finding method (secant, bisection, etc.)
23 vues (au cours des 30 derniers jours)
Afficher commentaires plus anciens
ABDALLA AL KHALEDI
le 18 Sep 2023
Commenté : ABDALLA AL KHALEDI
le 19 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
0 commentaires
Réponse acceptée
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);
1 commentaire
ABDALLA AL KHALEDI
le 19 Sep 2023
Modifié(e) : ABDALLA AL KHALEDI
le 19 Sep 2023
Plus de réponses (1)
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)
fprintf('Root k_hi = %.10f\n', k);
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
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.
Voir également
Catégories
En savoir plus sur Calculus dans Help Center et File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!