Programming wave dispersion equation

Hi everyone,

I am quite new to MATLAB and asking for some help to programme the wave dispersion equation.

ω^2= gktanh(kh)

where

𝑔 is gravity (9.81), ℎ is water depth, and 𝑘 is the wavenumber.

I would like to programme this function using an iterative process, and use it whenever I need to solve it for finding wavelength L, wave period T (2pi/omega)

Could someone please guide me in the process to create it and help further understand the mechanisms of iteration in MATLAB?

Thanks a lot!

[Torsten]

So given omega, g and h, you want to solve the equation for k ? Where are L and T in the equation ?

[GIULIA]

Yes, solve the equation for k.

And then hopefully understand the process so that I can re-arrange to solve for L and T :)

[Sam Chak]

Hi @GIULIA

If the wave number and the angular frequency ,

then

.

This is an implicit equation because the wavelength λ cannot be defined as a direct function of the other variables.

Try using the solve() command from the Symbolic Math Toolbox, as demonstrated here, before resorting to a numerical approach to solve the problem iteratively.

Answer 1

[Torsten]

Code taken from

https://uk.mathworks.com/matlabcentral/fileexchange/156009-wavelength-calculation-based-on-dispersion-relationship

% Parameters
T = 10;
h = 0.2;
g = 9.81;
tolerance = 1e-6;

% Initial guess for L
L0 = T^2/g;

% Fixed point iteration to compute L
Lsol = calculateWavelength(g, h, T, tolerance, L0);

% Plot curve
f = @(L) L - g * T^2 / (2 * pi) * tanh(2 * pi * h ./ L);
L = 1:0.1:20;
hold on
plot(L,f(L))
plot(Lsol,f(Lsol),'o')
hold off
grid on

function L = calculateWavelength(g, h, T, tolerance, L0)

    % Iterate until the tolerance is met
    while (1)
        L = g * T^2 / (2 * pi) * tanh(2 * pi * h / L0);

        if (abs(L - L0) < tolerance)
            break;
        end

        L0 = L;
    end
    
end

[GIULIA]

Thanks Torsten!

If i'd have to re-arrange to find T (knowing k and h) would it then be:

f = @(T) 2*pi / sqrt(g*k*tanh(k*h))

tolerance = 1e-6;

function T = calculatePeriod(g, h, tolerance, T0)

% Iterate until the tolerance is met

while (1)

T = (2 * pi) / sqrt(g * k * tanh (k*h));

if (abs(T - T0) < tolerance)

break;

end

T0 = T;

end

end

My questions in that case are:

1) What would I use as initial guess for T (T0)?

2) Can I always use 1e-6 as tolerance?

[Torsten]

You can explicitly solve for T:

T = sqrt( 2*pi*L/g * coth(2*pi*h/L) )

Answer 2

[Sam Chak]

An approach that directly use the solve() function.

syms lambda positive

% parameters
T   = 10;
h   = 0.2;
g   = 9.81;

% wave dispersion equation
eq  = lambda - g*T^2/(2*pi)*tanh(2*pi*h/lambda) == 0

eq = 

sol = solve(eq, lambda)

Warning: Unable to solve symbolically. Returning a numeric solution using vpasolve.

sol = 

13.988348869027691633555641776592

[GIULIA]

Thanks Sam!

I have tried also a similar approach using fsolve

function k= dispersion(T,h);

g= 9.81;

omega= 2*pi/ T;

%initialize k value

k = omega^2 / g;

dispf = @(k) omega^2 - g * k .* tanh(k * h);

no = optimset('Display','off'); %don't show output

k = fsolve(dispf, k, no);

My question is how I would re-arrange it to find for wave period t?

[Sam Chak]

Hi @GIULIA

Can you perform algebraic operations to rearrange the lambda (λ) equation to isolate the variable T on the left-hand side? This will allow you to verify whether you arrive at the same result as @Torsten in this comment.

[GIULIA]

Right, thank you!