I have the following code for my problem where i apply the match eigen function technique to find the unknown coefficients and then plotting one result wave motion inside a moonpool structure. However the expected plot for me should be decreasing after the resonance peak. But in my code it is decrease but then again blows up. I did not find any source of error. The analytical expressions which i used for my problem is write in pdf file at the end. Here is my code . I just want to figure out my source of error to fix it but i couldn’t. So I need help in acheiving my desired results. I also attached the refrence plot.

clear all;
close all;
clc;
tic;
% Parameters from your setup
N = 5; % Number of evanescent modes
M = 3; % Number of azimuthal modes
RE = 1.0; % Outer radius (m)
ratio = 2.0; % Ratio R_E / R_I
RI = RE / ratio; % Inner radius (m)
h = 4.0 * RE; % Water depth (m)
b = 2.0; % Distance from cylinder bottom to seabed (m)
d = h – b; % Interface depth (m)
g = 9.81; % Gravity (m/s^2)
l = 0; % Azimuthal order
% I_given_wavenumber = 1;
X_kRE = linspace(0.01, 4, 200);
% Initialize eta0 array before loop
eta0 = zeros(size(X_kRE)); % Preallocate
for idx = 1:length(X_kRE)
% if I_given_wavenumber == 1
wk = X_kRE(idx) / RE; % dimensional wavenumber
omega = sqrt(g * wk * tanh(wk * h)); % dispersion relation
fun_alpha = @(x) omega^2 + g*x.*tan(x*h); %dispersion equation for evanescent modes

for n = 1:N
if n == 1
k(n) = -1i * wk;
else

x0_left = (2*n – 3) * pi / (2*h);
x0_right = (2*n – 1) * pi / (2*h);
guess = (x0_left + x0_right)/2;

% Use lsqnonlin for better convergence
k(n) = lsqnonlin(fun_alpha, guess, x0_left, x0_right);
end
end
% % find evanescent modes k1…k_{N-1} stored in k_all(2)…k_all(N)
% for n = 2 : N
% m = n – 1; % mode index for evanescent
% a = ((2*m-1)*pi)/(2*h) + eps;
% b = (m*pi)/h – eps;
% % at a: tan→ -∞ so f(a)<0; at b: tan→0 so f(b)=+omega^2>0
% k_all(n) = fzero(@(k) g*k.*tan(k*h) + omega^2, [a, b]);
% end

% % % end
%% finding the root lemda_m =m*pi/b%%%%%%
for j=1:M
m(j)=(pi*(j-1))/b;

end
% %%%%%%%%%%%%%%%%%%%%% Define matrix A %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
A = zeros(M,N);
A(1,1) = (cosh(wk*h)*sinh(wk*b))/(wk*b*(2*wk*h+sinh(2*wk*h)))*(besselh(l, wk*RE)) /(dbesselh(l, wk*RE)); % Set A_00
for j = 2:N
A(1,j) =(cos(k(j)*h)*sin(k(j)*b) )/ (k(j)*b*(2*k(j)*h + sin(2*k(j)*h)))*( besselk(l, k(j)*RE) )…
/(dbesselk(l, k(j)*RE)); % Set A_0n

end

for i = 2:M
A(i,1) =(2 * cosh(wk*h) * (-1)^(i-1) * wk*sinh(wk*b)) / (b * (2*wk*h + sinh(2*wk*h))*(wk^2 + m(i)^2))…
* (besselh(l, wk*RE))/ (dbesselh(l, wk*RE)); % Set A_m0
end
for i = 2:M
for j = 2:N
A(i,j) = (2*cos(k(j)*h)*(-1)^(i-1) *k(j)* sin(k(j)*b))/ (b * (2*k(j)*h + sin(2*k(j)*h))*(k(j)^2-m(i)^2))…
*(besselk(l, k(j)*RE)) / (dbesselk(l, k(j)*RE));% Set A_mn
end
end
%%%%%%%%%%%%%%%%%%%%%%%% DefineMatrix B %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
B=zeros(N,M);
B(1,1) = (4*sinh(wk*b)) / (RE * wk*log(RE/RI) *cosh(wk*h)); %set B_00
%%%%%%%%%%%B_0m%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for j = 2:M
Rml_prime_RE = m(j)*(besselk(l, m(j)*RI) * dbesseli(l, m(j)*RE) …
– besseli(l, m(j)*RI) * dbesselk(l, m(j)*RE))…
/(besselk(l, m(j)*RI) * besseli(l, m(j)*RE) – besselk(l, m(j)*RE) * besseli(l, m(j)*RI));
B(1,j) = Rml_prime_RE * (4 * wk * (-1)^(j-1) * sinh(wk*b)) / (cosh(wk*h) * (wk^2 + m(j)^2));
end
%%%%%%%%%%%B_n0%%%%%%%%%%%%%%%%%%%%%%%%%%%
for i=2:N
B(i,1)=(4*sin(k(i)*b))/(RE*k(i)*log(RE/RI)*cos(k(i)*h));% Set B_n0
end
%%%%%%%%%%%%%%%%%%%%%%%%%%B_nm%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for i=2:N
for j=2:M
Rml_prime_RE = m(j)*(besselk(l, m(j)*RI) * dbesseli(l, m(j)*RE) …
– besseli(l, m(j)*RI) * dbesselk(l, m(j)*RE))…
/(besselk(l, m(j)*RI) * besseli(l, m(j)*RE) – besselk(l, m(j)*RE) * besseli(l, m(j)*RI));
B(i,j) = Rml_prime_RE * (4 * k(i) * (-1)^(j-1) * sin(k(i)*b)) / (cos(k(i)*h) * (k(i)^2 – m(j)^2));
end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Define Matrix C %%%%%%%%%%%%%%%%%%%%%%%%
C = zeros(N,M);
C(1,1) = -4 * sinh(wk*b) / (RE * wk*log(RE/RI) * cosh(wk*h));
% %%%%% DEfine C0m%%%%%%
for j = 2:M
Rml_prime_star_RE = m(j)*(besseli(l, m(j)*RE) * dbesselk(l, m(j)*RE) …
– besselk(l, m(j)*RE) * dbesseli(l, m(j)*RE))…
/(besselk(l, m(j)*RI) * besseli(l, m(j)*RE) – besselk(l, m(j)*RE) * besseli(l, m(j)*RI));
C(1,j) = Rml_prime_star_RE * (4 * wk * (-1)^(j-1) * sinh(wk*b)) / (cosh(wk*h) * (wk^2 + m(j)^2));
end
% %%%%%%% Define Cn0%%%%%%
for i = 2:N
C(i,1) = -4 * sin(k(i)*b) / (RE *k(i)* log(RE/RI) * cos(k(i)*h));
end
% % %%%%%% Define Cnm%%%%%%
for i = 2:N
for j =2:M
Rml_prime_star_RE = m(j)*(besseli(l, m(j)*RE) * dbesselk(l, m(j)*RE) …
– besselk(l, m(j)*RE) * dbesseli(l, m(j)*RE))…
/(besselk(l, m(j)*RI) * besseli(l, m(j)*RE) – besselk(l, m(j)*RE) * besseli(l, m(j)*RI));
C(i,j) = Rml_prime_star_RE * (4 * k(i) * (-1)^(j-1) * sin(k(i)*b)) / (cos(k(i)*h) * (k(i)^2 – m(j)^2));
end
end
%%%%%%% write Matrix D %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
D = zeros(M,N);
%%% write D00 %%%%%%
D(1,1) = (cosh(wk*h) * sinh(wk*b)) / (wk*b * (2*wk*h + sinh(2*wk*h))) * (besselj(l, wk*RI))/ (dbesselj(l, wk*RI));
%%%%% write D0n %%%%%%%%%%
for j =2:N
D(1,j) = (cos(k(j)*h) * sin(k(j)*b)) / (k(j)*b * (2*k(j)*h + sin(2*k(j)*h))) * (besseli(l, k(j)*RI) )/ (dbesseli(l, k(j)*RI));

end
%%%%%% write Dm0 %%%%%%%%%
for i = 2:M
D(i,1) = (2 * cosh(wk*h) * (-1)^(i-1) * wk * sinh(wk*b)) /(b * (2*wk*h + sinh(2*wk*h)) * (wk^2 + m(i)^2))…
*(besselj(l, wk*RI) )/(dbesselj(l, wk*RI));
end
%%%%% Define Dmn%%%%%%%%
for i = 2:M
for j = 2:N
D(i,j) = (2 * cos(k(j)*h) * (-1)^(i-1) * k(j) * sin(k(j)*b)) /(b * (2*k(j)*h + sin(2*k(j)*h)) * (k(j)^2 – m(i)^2))…
*(besseli(l, k(j)*RI)) / (dbesseli(l, k(j)*RI));
end
end
%%%%%%%%%%%%%%%%%%%%%%%%%% Define Matrix E %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
E = zeros(N, M); % Preallocate the matrix E
%%%%% Define E00%%%%%%%%%%%%
E(1,1) = (4 * sinh(wk*b)) / (RI *wk* log(RE/RI) * cosh(wk*h));
%%%% Define Eom %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for j = 2:M
Rml_prime_RI = m(j)*(besselk(l, m(j)*RI) * dbesseli(l, m(j)*RI) …
– besseli(l, m(j)*RI) * dbesselk(l, m(j)*RI))…
/ (besselk(l, m(j)*RI) * besseli(l, m(j)*RE) – besselk(l, m(j)*RE) * besseli(l, m(j)*RI));
E(1,j) = Rml_prime_RI * (4 * wk * (-1)^(j-1) * sinh(wk*b)) / (cosh(wk*h) * (wk^2 + m(j)^2));
%
end
%%%%%% Define Eno%%%%%%%%%%%%%%
for i = 2:N
E(i,1) = (4 * sin(k(i)*b)) / (RI *k(i)* log(RE/RI) * cos(k(i)*h));
end
for i = 2:N
for j = 2:M
Rml_prime_RI = m(j)*(besselk(l, m(j)*RI) * dbesseli(l, m(j)*RI) …
– besseli(l, m(j)*RI) * dbesselk(l, m(j)*RI))…
/ (besselk(l, m(j)*RI) * besseli(l, m(j)*RE) – besselk(l, m(j)*RE) * besseli(l, m(j)*RI));
E(i,j) = Rml_prime_RI * (4 * k(i) * (-1)^(j-1) * sin(k(i)*b)) / (cos(k(i)*h) * (k(i)^2 – m(j)^2));
end
end
%%%%%%%%%%%%%%%%%%%%%%%% Now write matrix F %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
F = zeros(N,M);
%%%%%%%% Define F00 %%%%%%%%%%
F(1,1) = (-4 * sinh(wk*b)) / (RI * wk*log(RE/RI) * cosh(wk*h));
% %%%%%% Define F0m%%%%
for j = 2:M
Rml_star_prime_RI = m(j)*(besseli(l, m(j)*RE) * dbesselk(l, m(j)*RI) …
– besselk(l, m(j)*RE) * dbesseli(l, m(j)*RI))/((besselk(l, m(j)*RI) …
* besseli(l, m(j)*RE) – besselk(l, m(j)*RE) * besseli(l, m(j)*RI)));
F(1,j) = Rml_star_prime_RI * (4 * wk * (-1)^(j-1) * sinh(wk*b)) / (cosh(wk*h) * (wk^2 + m(j)^2));
end
%%%%% Defin Fn0%%%%%%
for i = 2:N
F(i,1) = (-4 * sin(k(i)*b)) / (RI *k(i)* log(RE/RI) * cos(k(i)*h));
end
%%%%%%% Define Fnm %%%%%%%
for i = 2:N
for j = 2:M
Rml_star_prime_RI = m(j)*(besseli(l, m(j)*RE) * dbesselk(l, m(j)*RI) …
– besselk(l, m(j)*RE) * dbesseli(l, m(j)*RI))/((besselk(l, m(j)*RI) …
* besseli(l, m(j)*RE) – besselk(l, m(j)*RE) * besseli(l, m(j)*RI)));
F(i,j) = Rml_star_prime_RI * (4 * k(i)* (-1)^(j-1) * sin(k(i)*b)) / (cos(k(i)*h) * (k(i)^2 – m(j)^2));
end
end
% % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% % Define right-hand side vector G

G = zeros(2*M+2*N, 1); % Total length = M + N+ M+ N = 2M+2N
%
% Block 1: G (size M = 4)
for i = 1:M
if i == 1
G(i, 1) = (sinh(wk*b))/((b*wk)*cosh(wk*h))*besselj(l, wk*RE) ; % Z_0^l
else
G(i, 1) = (2 *wk*(-1)^(i-1)*sinh(wk*b))/(b *(wk^2 +m(i)^2))*besselj(l, wk*RE); %Z_m^l
end
end
%
% Block 2: Y (size N = 3)
for j = 1:N
if j == 1
G(j + M, 1) = -dbesselj(l, wk*RE) * (2*wk*h + sinh(2*wk*h)) /(cosh(wk*h)^2); % Y_0^l
else
G(j + M, 1) = 0; % Y_n^l
end
end

% Block 3: X (size M = 4)
for i = 1:M
G(i + M + N, 1) = 0; % X_0^l, X_m^l
end

% Block 4: W (size N = 3)
for j = 1:N
G(j + M + N + M, 1) = 0; % W_0^l, W_n^l
end
% Identity matrices
I_M = eye(M);
I_N = eye(N);

% Zero matrices
ZMM = zeros(M, M);
ZMN = zeros(M, N);
ZNM = zeros(N, M);
ZNN = zeros(N, N);

% Construct the full block matrix H
H = [ I_M, -A, ZMM, ZMN; % Block row 1
-B, I_N, -C, ZNN; % Block row 2
ZMM, ZMN, I_M, -D; % Block row 3
-E, ZNN, -F, I_N]; % Block row 4
% X = H G; % Solve the linear system
X = pinv(H) * G; % Use pseudoinverse for stability
%
b_vec = X(1:M); % Coefficients b^l
a_vec = X(M+1 : M+N); % Coefficients a^l
c_vec = X(M+N+1 : 2*M+N); % Coefficients c^l
d_vec = X(2*M+N+1 : end); % Coefficients d^l
%%%%%%%%% Wave motion%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
term1 = (d_vec(1)*cosh(wk*h)^2) / ((2*wk*h + sinh(2*wk*h)) * dbesselj(l, wk*RI));
sum1 = 0;
for n =2:N
sum1 = sum1 + d_vec(n)*cos(k(n)*h)^2/ ((2*k(n)*h + sin(2*k(n)*h))* dbesseli(l, k(n)*RI));
end
eta0(idx) = abs(term1+sum1);
% % disp(k.’) % Should all be real and positive for n ≥ 2
% %
% %

end

plot(X_kRE, eta0, ‘k’, ‘LineWidth’, 1.5);
xlabel(‘$k_0 R_E$’, ‘Interpreter’, ‘latex’);
ylabel(‘$eta / (iA)$’, ‘Interpreter’, ‘latex’);
title(‘Wave motion amplitude for RE/RI = 4.0’);
grid on;
xlim([0 4]); % Match the reference plot
ylim([0 6]); % Optional: based on expected peak height
elapsedTime = toc;
disp([‘Time consuming = ‘, num2str(elapsedTime), ‘ s’]);

function out = dbesselh(l, z)
if l == 0
out = -besselh(1, 1, z);
else
out = 0.5 * (besselh(l – 1, 1, z) – besselh(l + 1, 1, z));
end
end

function out = dbesseli(l, z)
if l == 0
out = besseli(1, z);
else
out = 0.5 * (besseli(l – 1, z) + besseli(l + 1, z));
end
end

function out = dbesselj(l, z)
if l == 0
out = -besselj(1, z);
else
out = 0.5 * (besselj(l – 1, z) – besselj(l + 1, z));
end
end

function out = dbesselk(l, z)
if l == 0
out = -besselk(1, z);
else
out = -0.5 * (besselk(l – 1, z) + besselk(l + 1, z));
end
endI have the following code for my problem where i apply the match eigen function technique to find the unknown coefficients and then plotting one result wave motion inside a moonpool structure. However the expected plot for me should be decreasing after the resonance peak. But in my code it is decrease but then again blows up. I did not find any source of error. The analytical expressions which i used for my problem is write in pdf file at the end. Here is my code . I just want to figure out my source of error to fix it but i couldn’t. So I need help in acheiving my desired results. I also attached the refrence plot.

clear all;
close all;
clc;
tic;
% Parameters from your setup
N = 5; % Number of evanescent modes
M = 3; % Number of azimuthal modes
RE = 1.0; % Outer radius (m)
ratio = 2.0; % Ratio R_E / R_I
RI = RE / ratio; % Inner radius (m)
h = 4.0 * RE; % Water depth (m)
b = 2.0; % Distance from cylinder bottom to seabed (m)
d = h – b; % Interface depth (m)
g = 9.81; % Gravity (m/s^2)
l = 0; % Azimuthal order
% I_given_wavenumber = 1;
X_kRE = linspace(0.01, 4, 200);
% Initialize eta0 array before loop
eta0 = zeros(size(X_kRE)); % Preallocate
for idx = 1:length(X_kRE)
% if I_given_wavenumber == 1
wk = X_kRE(idx) / RE; % dimensional wavenumber
omega = sqrt(g * wk * tanh(wk * h)); % dispersion relation
fun_alpha = @(x) omega^2 + g*x.*tan(x*h); %dispersion equation for evanescent modes

for n = 1:N
if n == 1
k(n) = -1i * wk;
else

x0_left = (2*n – 3) * pi / (2*h);
x0_right = (2*n – 1) * pi / (2*h);
guess = (x0_left + x0_right)/2;

% Use lsqnonlin for better convergence
k(n) = lsqnonlin(fun_alpha, guess, x0_left, x0_right);
end
end
% % find evanescent modes k1…k_{N-1} stored in k_all(2)…k_all(N)
% for n = 2 : N
% m = n – 1; % mode index for evanescent
% a = ((2*m-1)*pi)/(2*h) + eps;
% b = (m*pi)/h – eps;
% % at a: tan→ -∞ so f(a)<0; at b: tan→0 so f(b)=+omega^2>0
% k_all(n) = fzero(@(k) g*k.*tan(k*h) + omega^2, [a, b]);
% end

% % % end
%% finding the root lemda_m =m*pi/b%%%%%%
for j=1:M
m(j)=(pi*(j-1))/b;

end
% %%%%%%%%%%%%%%%%%%%%% Define matrix A %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
A = zeros(M,N);
A(1,1) = (cosh(wk*h)*sinh(wk*b))/(wk*b*(2*wk*h+sinh(2*wk*h)))*(besselh(l, wk*RE)) /(dbesselh(l, wk*RE)); % Set A_00
for j = 2:N
A(1,j) =(cos(k(j)*h)*sin(k(j)*b) )/ (k(j)*b*(2*k(j)*h + sin(2*k(j)*h)))*( besselk(l, k(j)*RE) )…
/(dbesselk(l, k(j)*RE)); % Set A_0n

end

for i = 2:M
A(i,1) =(2 * cosh(wk*h) * (-1)^(i-1) * wk*sinh(wk*b)) / (b * (2*wk*h + sinh(2*wk*h))*(wk^2 + m(i)^2))…
* (besselh(l, wk*RE))/ (dbesselh(l, wk*RE)); % Set A_m0
end
for i = 2:M
for j = 2:N
A(i,j) = (2*cos(k(j)*h)*(-1)^(i-1) *k(j)* sin(k(j)*b))/ (b * (2*k(j)*h + sin(2*k(j)*h))*(k(j)^2-m(i)^2))…
*(besselk(l, k(j)*RE)) / (dbesselk(l, k(j)*RE));% Set A_mn
end
end
%%%%%%%%%%%%%%%%%%%%%%%% DefineMatrix B %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
B=zeros(N,M);
B(1,1) = (4*sinh(wk*b)) / (RE * wk*log(RE/RI) *cosh(wk*h)); %set B_00
%%%%%%%%%%%B_0m%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for j = 2:M
Rml_prime_RE = m(j)*(besselk(l, m(j)*RI) * dbesseli(l, m(j)*RE) …
– besseli(l, m(j)*RI) * dbesselk(l, m(j)*RE))…
/(besselk(l, m(j)*RI) * besseli(l, m(j)*RE) – besselk(l, m(j)*RE) * besseli(l, m(j)*RI));
B(1,j) = Rml_prime_RE * (4 * wk * (-1)^(j-1) * sinh(wk*b)) / (cosh(wk*h) * (wk^2 + m(j)^2));
end
%%%%%%%%%%%B_n0%%%%%%%%%%%%%%%%%%%%%%%%%%%
for i=2:N
B(i,1)=(4*sin(k(i)*b))/(RE*k(i)*log(RE/RI)*cos(k(i)*h));% Set B_n0
end
%%%%%%%%%%%%%%%%%%%%%%%%%%B_nm%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for i=2:N
for j=2:M
Rml_prime_RE = m(j)*(besselk(l, m(j)*RI) * dbesseli(l, m(j)*RE) …
– besseli(l, m(j)*RI) * dbesselk(l, m(j)*RE))…
/(besselk(l, m(j)*RI) * besseli(l, m(j)*RE) – besselk(l, m(j)*RE) * besseli(l, m(j)*RI));
B(i,j) = Rml_prime_RE * (4 * k(i) * (-1)^(j-1) * sin(k(i)*b)) / (cos(k(i)*h) * (k(i)^2 – m(j)^2));
end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Define Matrix C %%%%%%%%%%%%%%%%%%%%%%%%
C = zeros(N,M);
C(1,1) = -4 * sinh(wk*b) / (RE * wk*log(RE/RI) * cosh(wk*h));
% %%%%% DEfine C0m%%%%%%
for j = 2:M
Rml_prime_star_RE = m(j)*(besseli(l, m(j)*RE) * dbesselk(l, m(j)*RE) …
– besselk(l, m(j)*RE) * dbesseli(l, m(j)*RE))…
/(besselk(l, m(j)*RI) * besseli(l, m(j)*RE) – besselk(l, m(j)*RE) * besseli(l, m(j)*RI));
C(1,j) = Rml_prime_star_RE * (4 * wk * (-1)^(j-1) * sinh(wk*b)) / (cosh(wk*h) * (wk^2 + m(j)^2));
end
% %%%%%%% Define Cn0%%%%%%
for i = 2:N
C(i,1) = -4 * sin(k(i)*b) / (RE *k(i)* log(RE/RI) * cos(k(i)*h));
end
% % %%%%%% Define Cnm%%%%%%
for i = 2:N
for j =2:M
Rml_prime_star_RE = m(j)*(besseli(l, m(j)*RE) * dbesselk(l, m(j)*RE) …
– besselk(l, m(j)*RE) * dbesseli(l, m(j)*RE))…
/(besselk(l, m(j)*RI) * besseli(l, m(j)*RE) – besselk(l, m(j)*RE) * besseli(l, m(j)*RI));
C(i,j) = Rml_prime_star_RE * (4 * k(i) * (-1)^(j-1) * sin(k(i)*b)) / (cos(k(i)*h) * (k(i)^2 – m(j)^2));
end
end
%%%%%%% write Matrix D %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
D = zeros(M,N);
%%% write D00 %%%%%%
D(1,1) = (cosh(wk*h) * sinh(wk*b)) / (wk*b * (2*wk*h + sinh(2*wk*h))) * (besselj(l, wk*RI))/ (dbesselj(l, wk*RI));
%%%%% write D0n %%%%%%%%%%
for j =2:N
D(1,j) = (cos(k(j)*h) * sin(k(j)*b)) / (k(j)*b * (2*k(j)*h + sin(2*k(j)*h))) * (besseli(l, k(j)*RI) )/ (dbesseli(l, k(j)*RI));

end
%%%%%% write Dm0 %%%%%%%%%
for i = 2:M
D(i,1) = (2 * cosh(wk*h) * (-1)^(i-1) * wk * sinh(wk*b)) /(b * (2*wk*h + sinh(2*wk*h)) * (wk^2 + m(i)^2))…
*(besselj(l, wk*RI) )/(dbesselj(l, wk*RI));
end
%%%%% Define Dmn%%%%%%%%
for i = 2:M
for j = 2:N
D(i,j) = (2 * cos(k(j)*h) * (-1)^(i-1) * k(j) * sin(k(j)*b)) /(b * (2*k(j)*h + sin(2*k(j)*h)) * (k(j)^2 – m(i)^2))…
*(besseli(l, k(j)*RI)) / (dbesseli(l, k(j)*RI));
end
end
%%%%%%%%%%%%%%%%%%%%%%%%%% Define Matrix E %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
E = zeros(N, M); % Preallocate the matrix E
%%%%% Define E00%%%%%%%%%%%%
E(1,1) = (4 * sinh(wk*b)) / (RI *wk* log(RE/RI) * cosh(wk*h));
%%%% Define Eom %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for j = 2:M
Rml_prime_RI = m(j)*(besselk(l, m(j)*RI) * dbesseli(l, m(j)*RI) …
– besseli(l, m(j)*RI) * dbesselk(l, m(j)*RI))…
/ (besselk(l, m(j)*RI) * besseli(l, m(j)*RE) – besselk(l, m(j)*RE) * besseli(l, m(j)*RI));
E(1,j) = Rml_prime_RI * (4 * wk * (-1)^(j-1) * sinh(wk*b)) / (cosh(wk*h) * (wk^2 + m(j)^2));
%
end
%%%%%% Define Eno%%%%%%%%%%%%%%
for i = 2:N
E(i,1) = (4 * sin(k(i)*b)) / (RI *k(i)* log(RE/RI) * cos(k(i)*h));
end
for i = 2:N
for j = 2:M
Rml_prime_RI = m(j)*(besselk(l, m(j)*RI) * dbesseli(l, m(j)*RI) …
– besseli(l, m(j)*RI) * dbesselk(l, m(j)*RI))…
/ (besselk(l, m(j)*RI) * besseli(l, m(j)*RE) – besselk(l, m(j)*RE) * besseli(l, m(j)*RI));
E(i,j) = Rml_prime_RI * (4 * k(i) * (-1)^(j-1) * sin(k(i)*b)) / (cos(k(i)*h) * (k(i)^2 – m(j)^2));
end
end
%%%%%%%%%%%%%%%%%%%%%%%% Now write matrix F %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
F = zeros(N,M);
%%%%%%%% Define F00 %%%%%%%%%%
F(1,1) = (-4 * sinh(wk*b)) / (RI * wk*log(RE/RI) * cosh(wk*h));
% %%%%%% Define F0m%%%%
for j = 2:M
Rml_star_prime_RI = m(j)*(besseli(l, m(j)*RE) * dbesselk(l, m(j)*RI) …
– besselk(l, m(j)*RE) * dbesseli(l, m(j)*RI))/((besselk(l, m(j)*RI) …
* besseli(l, m(j)*RE) – besselk(l, m(j)*RE) * besseli(l, m(j)*RI)));
F(1,j) = Rml_star_prime_RI * (4 * wk * (-1)^(j-1) * sinh(wk*b)) / (cosh(wk*h) * (wk^2 + m(j)^2));
end
%%%%% Defin Fn0%%%%%%
for i = 2:N
F(i,1) = (-4 * sin(k(i)*b)) / (RI *k(i)* log(RE/RI) * cos(k(i)*h));
end
%%%%%%% Define Fnm %%%%%%%
for i = 2:N
for j = 2:M
Rml_star_prime_RI = m(j)*(besseli(l, m(j)*RE) * dbesselk(l, m(j)*RI) …
– besselk(l, m(j)*RE) * dbesseli(l, m(j)*RI))/((besselk(l, m(j)*RI) …
* besseli(l, m(j)*RE) – besselk(l, m(j)*RE) * besseli(l, m(j)*RI)));
F(i,j) = Rml_star_prime_RI * (4 * k(i)* (-1)^(j-1) * sin(k(i)*b)) / (cos(k(i)*h) * (k(i)^2 – m(j)^2));
end
end
% % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% % Define right-hand side vector G

G = zeros(2*M+2*N, 1); % Total length = M + N+ M+ N = 2M+2N
%
% Block 1: G (size M = 4)
for i = 1:M
if i == 1
G(i, 1) = (sinh(wk*b))/((b*wk)*cosh(wk*h))*besselj(l, wk*RE) ; % Z_0^l
else
G(i, 1) = (2 *wk*(-1)^(i-1)*sinh(wk*b))/(b *(wk^2 +m(i)^2))*besselj(l, wk*RE); %Z_m^l
end
end
%
% Block 2: Y (size N = 3)
for j = 1:N
if j == 1
G(j + M, 1) = -dbesselj(l, wk*RE) * (2*wk*h + sinh(2*wk*h)) /(cosh(wk*h)^2); % Y_0^l
else
G(j + M, 1) = 0; % Y_n^l
end
end

% Block 3: X (size M = 4)
for i = 1:M
G(i + M + N, 1) = 0; % X_0^l, X_m^l
end

% Block 4: W (size N = 3)
for j = 1:N
G(j + M + N + M, 1) = 0; % W_0^l, W_n^l
end
% Identity matrices
I_M = eye(M);
I_N = eye(N);

% Zero matrices
ZMM = zeros(M, M);
ZMN = zeros(M, N);
ZNM = zeros(N, M);
ZNN = zeros(N, N);

% Construct the full block matrix H
H = [ I_M, -A, ZMM, ZMN; % Block row 1
-B, I_N, -C, ZNN; % Block row 2
ZMM, ZMN, I_M, -D; % Block row 3
-E, ZNN, -F, I_N]; % Block row 4
% X = H G; % Solve the linear system
X = pinv(H) * G; % Use pseudoinverse for stability
%
b_vec = X(1:M); % Coefficients b^l
a_vec = X(M+1 : M+N); % Coefficients a^l
c_vec = X(M+N+1 : 2*M+N); % Coefficients c^l
d_vec = X(2*M+N+1 : end); % Coefficients d^l
%%%%%%%%% Wave motion%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
term1 = (d_vec(1)*cosh(wk*h)^2) / ((2*wk*h + sinh(2*wk*h)) * dbesselj(l, wk*RI));
sum1 = 0;
for n =2:N
sum1 = sum1 + d_vec(n)*cos(k(n)*h)^2/ ((2*k(n)*h + sin(2*k(n)*h))* dbesseli(l, k(n)*RI));
end
eta0(idx) = abs(term1+sum1);
% % disp(k.’) % Should all be real and positive for n ≥ 2
% %
% %

end

plot(X_kRE, eta0, ‘k’, ‘LineWidth’, 1.5);
xlabel(‘$k_0 R_E$’, ‘Interpreter’, ‘latex’);
ylabel(‘$eta / (iA)$’, ‘Interpreter’, ‘latex’);
title(‘Wave motion amplitude for RE/RI = 4.0’);
grid on;
xlim([0 4]); % Match the reference plot
ylim([0 6]); % Optional: based on expected peak height
elapsedTime = toc;
disp([‘Time consuming = ‘, num2str(elapsedTime), ‘ s’]);

function out = dbesselh(l, z)
if l == 0
out = -besselh(1, 1, z);
else
out = 0.5 * (besselh(l – 1, 1, z) – besselh(l + 1, 1, z));
end
end

function out = dbesseli(l, z)
if l == 0
out = besseli(1, z);
else
out = 0.5 * (besseli(l – 1, z) + besseli(l + 1, z));
end
end

function out = dbesselj(l, z)
if l == 0
out = -besselj(1, z);
else
out = 0.5 * (besselj(l – 1, z) – besselj(l + 1, z));
end
end

function out = dbesselk(l, z)
if l == 0
out = -besselk(1, z);
else
out = -0.5 * (besselk(l – 1, z) + besselk(l + 1, z));
end
end I have the following code for my problem where i apply the match eigen function technique to find the unknown coefficients and then plotting one result wave motion inside a moonpool structure. However the expected plot for me should be decreasing after the resonance peak. But in my code it is decrease but then again blows up. I did not find any source of error. The analytical expressions which i used for my problem is write in pdf file at the end. Here is my code . I just want to figure out my source of error to fix it but i couldn’t. So I need help in acheiving my desired results. I also attached the refrence plot.

clear all;
close all;
clc;
tic;
% Parameters from your setup
N = 5; % Number of evanescent modes
M = 3; % Number of azimuthal modes
RE = 1.0; % Outer radius (m)
ratio = 2.0; % Ratio R_E / R_I
RI = RE / ratio; % Inner radius (m)
h = 4.0 * RE; % Water depth (m)
b = 2.0; % Distance from cylinder bottom to seabed (m)
d = h – b; % Interface depth (m)
g = 9.81; % Gravity (m/s^2)
l = 0; % Azimuthal order
% I_given_wavenumber = 1;
X_kRE = linspace(0.01, 4, 200);
% Initialize eta0 array before loop
eta0 = zeros(size(X_kRE)); % Preallocate
for idx = 1:length(X_kRE)
% if I_given_wavenumber == 1
wk = X_kRE(idx) / RE; % dimensional wavenumber
omega = sqrt(g * wk * tanh(wk * h)); % dispersion relation
fun_alpha = @(x) omega^2 + g*x.*tan(x*h); %dispersion equation for evanescent modes

for n = 1:N
if n == 1
k(n) = -1i * wk;
else

x0_left = (2*n – 3) * pi / (2*h);
x0_right = (2*n – 1) * pi / (2*h);
guess = (x0_left + x0_right)/2;

% Use lsqnonlin for better convergence
k(n) = lsqnonlin(fun_alpha, guess, x0_left, x0_right);
end
end
% % find evanescent modes k1…k_{N-1} stored in k_all(2)…k_all(N)
% for n = 2 : N
% m = n – 1; % mode index for evanescent
% a = ((2*m-1)*pi)/(2*h) + eps;
% b = (m*pi)/h – eps;
% % at a: tan→ -∞ so f(a)<0; at b: tan→0 so f(b)=+omega^2>0
% k_all(n) = fzero(@(k) g*k.*tan(k*h) + omega^2, [a, b]);
% end

% % % end
%% finding the root lemda_m =m*pi/b%%%%%%
for j=1:M
m(j)=(pi*(j-1))/b;

end
% %%%%%%%%%%%%%%%%%%%%% Define matrix A %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
A = zeros(M,N);
A(1,1) = (cosh(wk*h)*sinh(wk*b))/(wk*b*(2*wk*h+sinh(2*wk*h)))*(besselh(l, wk*RE)) /(dbesselh(l, wk*RE)); % Set A_00
for j = 2:N
A(1,j) =(cos(k(j)*h)*sin(k(j)*b) )/ (k(j)*b*(2*k(j)*h + sin(2*k(j)*h)))*( besselk(l, k(j)*RE) )…
/(dbesselk(l, k(j)*RE)); % Set A_0n

end

for i = 2:M
A(i,1) =(2 * cosh(wk*h) * (-1)^(i-1) * wk*sinh(wk*b)) / (b * (2*wk*h + sinh(2*wk*h))*(wk^2 + m(i)^2))…
* (besselh(l, wk*RE))/ (dbesselh(l, wk*RE)); % Set A_m0
end
for i = 2:M
for j = 2:N
A(i,j) = (2*cos(k(j)*h)*(-1)^(i-1) *k(j)* sin(k(j)*b))/ (b * (2*k(j)*h + sin(2*k(j)*h))*(k(j)^2-m(i)^2))…
*(besselk(l, k(j)*RE)) / (dbesselk(l, k(j)*RE));% Set A_mn
end
end
%%%%%%%%%%%%%%%%%%%%%%%% DefineMatrix B %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
B=zeros(N,M);
B(1,1) = (4*sinh(wk*b)) / (RE * wk*log(RE/RI) *cosh(wk*h)); %set B_00
%%%%%%%%%%%B_0m%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for j = 2:M
Rml_prime_RE = m(j)*(besselk(l, m(j)*RI) * dbesseli(l, m(j)*RE) …
– besseli(l, m(j)*RI) * dbesselk(l, m(j)*RE))…
/(besselk(l, m(j)*RI) * besseli(l, m(j)*RE) – besselk(l, m(j)*RE) * besseli(l, m(j)*RI));
B(1,j) = Rml_prime_RE * (4 * wk * (-1)^(j-1) * sinh(wk*b)) / (cosh(wk*h) * (wk^2 + m(j)^2));
end
%%%%%%%%%%%B_n0%%%%%%%%%%%%%%%%%%%%%%%%%%%
for i=2:N
B(i,1)=(4*sin(k(i)*b))/(RE*k(i)*log(RE/RI)*cos(k(i)*h));% Set B_n0
end
%%%%%%%%%%%%%%%%%%%%%%%%%%B_nm%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for i=2:N
for j=2:M
Rml_prime_RE = m(j)*(besselk(l, m(j)*RI) * dbesseli(l, m(j)*RE) …
– besseli(l, m(j)*RI) * dbesselk(l, m(j)*RE))…
/(besselk(l, m(j)*RI) * besseli(l, m(j)*RE) – besselk(l, m(j)*RE) * besseli(l, m(j)*RI));
B(i,j) = Rml_prime_RE * (4 * k(i) * (-1)^(j-1) * sin(k(i)*b)) / (cos(k(i)*h) * (k(i)^2 – m(j)^2));
end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Define Matrix C %%%%%%%%%%%%%%%%%%%%%%%%
C = zeros(N,M);
C(1,1) = -4 * sinh(wk*b) / (RE * wk*log(RE/RI) * cosh(wk*h));
% %%%%% DEfine C0m%%%%%%
for j = 2:M
Rml_prime_star_RE = m(j)*(besseli(l, m(j)*RE) * dbesselk(l, m(j)*RE) …
– besselk(l, m(j)*RE) * dbesseli(l, m(j)*RE))…
/(besselk(l, m(j)*RI) * besseli(l, m(j)*RE) – besselk(l, m(j)*RE) * besseli(l, m(j)*RI));
C(1,j) = Rml_prime_star_RE * (4 * wk * (-1)^(j-1) * sinh(wk*b)) / (cosh(wk*h) * (wk^2 + m(j)^2));
end
% %%%%%%% Define Cn0%%%%%%
for i = 2:N
C(i,1) = -4 * sin(k(i)*b) / (RE *k(i)* log(RE/RI) * cos(k(i)*h));
end
% % %%%%%% Define Cnm%%%%%%
for i = 2:N
for j =2:M
Rml_prime_star_RE = m(j)*(besseli(l, m(j)*RE) * dbesselk(l, m(j)*RE) …
– besselk(l, m(j)*RE) * dbesseli(l, m(j)*RE))…
/(besselk(l, m(j)*RI) * besseli(l, m(j)*RE) – besselk(l, m(j)*RE) * besseli(l, m(j)*RI));
C(i,j) = Rml_prime_star_RE * (4 * k(i) * (-1)^(j-1) * sin(k(i)*b)) / (cos(k(i)*h) * (k(i)^2 – m(j)^2));
end
end
%%%%%%% write Matrix D %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
D = zeros(M,N);
%%% write D00 %%%%%%
D(1,1) = (cosh(wk*h) * sinh(wk*b)) / (wk*b * (2*wk*h + sinh(2*wk*h))) * (besselj(l, wk*RI))/ (dbesselj(l, wk*RI));
%%%%% write D0n %%%%%%%%%%
for j =2:N
D(1,j) = (cos(k(j)*h) * sin(k(j)*b)) / (k(j)*b * (2*k(j)*h + sin(2*k(j)*h))) * (besseli(l, k(j)*RI) )/ (dbesseli(l, k(j)*RI));

end
%%%%%% write Dm0 %%%%%%%%%
for i = 2:M
D(i,1) = (2 * cosh(wk*h) * (-1)^(i-1) * wk * sinh(wk*b)) /(b * (2*wk*h + sinh(2*wk*h)) * (wk^2 + m(i)^2))…
*(besselj(l, wk*RI) )/(dbesselj(l, wk*RI));
end
%%%%% Define Dmn%%%%%%%%
for i = 2:M
for j = 2:N
D(i,j) = (2 * cos(k(j)*h) * (-1)^(i-1) * k(j) * sin(k(j)*b)) /(b * (2*k(j)*h + sin(2*k(j)*h)) * (k(j)^2 – m(i)^2))…
*(besseli(l, k(j)*RI)) / (dbesseli(l, k(j)*RI));
end
end
%%%%%%%%%%%%%%%%%%%%%%%%%% Define Matrix E %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
E = zeros(N, M); % Preallocate the matrix E
%%%%% Define E00%%%%%%%%%%%%
E(1,1) = (4 * sinh(wk*b)) / (RI *wk* log(RE/RI) * cosh(wk*h));
%%%% Define Eom %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for j = 2:M
Rml_prime_RI = m(j)*(besselk(l, m(j)*RI) * dbesseli(l, m(j)*RI) …
– besseli(l, m(j)*RI) * dbesselk(l, m(j)*RI))…
/ (besselk(l, m(j)*RI) * besseli(l, m(j)*RE) – besselk(l, m(j)*RE) * besseli(l, m(j)*RI));
E(1,j) = Rml_prime_RI * (4 * wk * (-1)^(j-1) * sinh(wk*b)) / (cosh(wk*h) * (wk^2 + m(j)^2));
%
end
%%%%%% Define Eno%%%%%%%%%%%%%%
for i = 2:N
E(i,1) = (4 * sin(k(i)*b)) / (RI *k(i)* log(RE/RI) * cos(k(i)*h));
end
for i = 2:N
for j = 2:M
Rml_prime_RI = m(j)*(besselk(l, m(j)*RI) * dbesseli(l, m(j)*RI) …
– besseli(l, m(j)*RI) * dbesselk(l, m(j)*RI))…
/ (besselk(l, m(j)*RI) * besseli(l, m(j)*RE) – besselk(l, m(j)*RE) * besseli(l, m(j)*RI));
E(i,j) = Rml_prime_RI * (4 * k(i) * (-1)^(j-1) * sin(k(i)*b)) / (cos(k(i)*h) * (k(i)^2 – m(j)^2));
end
end
%%%%%%%%%%%%%%%%%%%%%%%% Now write matrix F %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
F = zeros(N,M);
%%%%%%%% Define F00 %%%%%%%%%%
F(1,1) = (-4 * sinh(wk*b)) / (RI * wk*log(RE/RI) * cosh(wk*h));
% %%%%%% Define F0m%%%%
for j = 2:M
Rml_star_prime_RI = m(j)*(besseli(l, m(j)*RE) * dbesselk(l, m(j)*RI) …
– besselk(l, m(j)*RE) * dbesseli(l, m(j)*RI))/((besselk(l, m(j)*RI) …
* besseli(l, m(j)*RE) – besselk(l, m(j)*RE) * besseli(l, m(j)*RI)));
F(1,j) = Rml_star_prime_RI * (4 * wk * (-1)^(j-1) * sinh(wk*b)) / (cosh(wk*h) * (wk^2 + m(j)^2));
end
%%%%% Defin Fn0%%%%%%
for i = 2:N
F(i,1) = (-4 * sin(k(i)*b)) / (RI *k(i)* log(RE/RI) * cos(k(i)*h));
end
%%%%%%% Define Fnm %%%%%%%
for i = 2:N
for j = 2:M
Rml_star_prime_RI = m(j)*(besseli(l, m(j)*RE) * dbesselk(l, m(j)*RI) …
– besselk(l, m(j)*RE) * dbesseli(l, m(j)*RI))/((besselk(l, m(j)*RI) …
* besseli(l, m(j)*RE) – besselk(l, m(j)*RE) * besseli(l, m(j)*RI)));
F(i,j) = Rml_star_prime_RI * (4 * k(i)* (-1)^(j-1) * sin(k(i)*b)) / (cos(k(i)*h) * (k(i)^2 – m(j)^2));
end
end
% % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% % Define right-hand side vector G

G = zeros(2*M+2*N, 1); % Total length = M + N+ M+ N = 2M+2N
%
% Block 1: G (size M = 4)
for i = 1:M
if i == 1
G(i, 1) = (sinh(wk*b))/((b*wk)*cosh(wk*h))*besselj(l, wk*RE) ; % Z_0^l
else
G(i, 1) = (2 *wk*(-1)^(i-1)*sinh(wk*b))/(b *(wk^2 +m(i)^2))*besselj(l, wk*RE); %Z_m^l
end
end
%
% Block 2: Y (size N = 3)
for j = 1:N
if j == 1
G(j + M, 1) = -dbesselj(l, wk*RE) * (2*wk*h + sinh(2*wk*h)) /(cosh(wk*h)^2); % Y_0^l
else
G(j + M, 1) = 0; % Y_n^l
end
end

% Block 3: X (size M = 4)
for i = 1:M
G(i + M + N, 1) = 0; % X_0^l, X_m^l
end

% Block 4: W (size N = 3)
for j = 1:N
G(j + M + N + M, 1) = 0; % W_0^l, W_n^l
end
% Identity matrices
I_M = eye(M);
I_N = eye(N);

% Zero matrices
ZMM = zeros(M, M);
ZMN = zeros(M, N);
ZNM = zeros(N, M);
ZNN = zeros(N, N);

% Construct the full block matrix H
H = [ I_M, -A, ZMM, ZMN; % Block row 1
-B, I_N, -C, ZNN; % Block row 2
ZMM, ZMN, I_M, -D; % Block row 3
-E, ZNN, -F, I_N]; % Block row 4
% X = H G; % Solve the linear system
X = pinv(H) * G; % Use pseudoinverse for stability
%
b_vec = X(1:M); % Coefficients b^l
a_vec = X(M+1 : M+N); % Coefficients a^l
c_vec = X(M+N+1 : 2*M+N); % Coefficients c^l
d_vec = X(2*M+N+1 : end); % Coefficients d^l
%%%%%%%%% Wave motion%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
term1 = (d_vec(1)*cosh(wk*h)^2) / ((2*wk*h + sinh(2*wk*h)) * dbesselj(l, wk*RI));
sum1 = 0;
for n =2:N
sum1 = sum1 + d_vec(n)*cos(k(n)*h)^2/ ((2*k(n)*h + sin(2*k(n)*h))* dbesseli(l, k(n)*RI));
end
eta0(idx) = abs(term1+sum1);
% % disp(k.’) % Should all be real and positive for n ≥ 2
% %
% %

end

plot(X_kRE, eta0, ‘k’, ‘LineWidth’, 1.5);
xlabel(‘$k_0 R_E$’, ‘Interpreter’, ‘latex’);
ylabel(‘$eta / (iA)$’, ‘Interpreter’, ‘latex’);
title(‘Wave motion amplitude for RE/RI = 4.0’);
grid on;
xlim([0 4]); % Match the reference plot
ylim([0 6]); % Optional: based on expected peak height
elapsedTime = toc;
disp([‘Time consuming = ‘, num2str(elapsedTime), ‘ s’]);

function out = dbesselh(l, z)
if l == 0
out = -besselh(1, 1, z);
else
out = 0.5 * (besselh(l – 1, 1, z) – besselh(l + 1, 1, z));
end
end

function out = dbesseli(l, z)
if l == 0
out = besseli(1, z);
else
out = 0.5 * (besseli(l – 1, z) + besseli(l + 1, z));
end
end

function out = dbesselj(l, z)
if l == 0
out = -besselj(1, z);
else
out = 0.5 * (besselj(l – 1, z) – besselj(l + 1, z));
end
end

function out = dbesselk(l, z)
if l == 0
out = -besselk(1, z);
else
out = -0.5 * (besselk(l – 1, z) + besselk(l + 1, z));
end
end moonpool, resonance, system of linear equation, semi analytical solution, wave motion MATLAB Answers — New Questions

​

I am working on a project where I have to analyse high frequency wind data (50 Hz) coming from a wind turbine data measurement system. this amount of data must be converted from .dat binary files to .mat files which I can use in matlab. The data must then be filtered and then averaged over 10 minutes to be compared to the data from another measurement system. Doing all this requires the analysis of thousands of data and it’s very time consuming (right now it is about 105 seconds just for the data of 2 days). How can I speed the process up?
for ii = 3:length(filedir)
filename = filedir(ii).name;
newdata = ReadFamosDataIntoTimeTable(filename);
% filter
IsConsidered = newdata.avbladeangleGRe<40 … % normal operation
& newdata.RAWS9>0.5 … % good availability
& ~isnan(newdata.iv10mswindspeed2GRe) & newdata.av100msabswinddirectionGRe>180 & newdata.av100msabswinddirectionGRe<250;
TurbineData(ii-2).date = filename;
TurbineData(ii-2).iv10mswindspeed2GRe = mean(newdata.iv10mswindspeed2GRe(IsConsidered));
TurbineData(ii-2).ivactivepowerGRe = mean(newdata.ivactivepowerGRe(IsConsidered));
TurbineData(ii-2).CalculatedAirdensity_GRe = mean(newdata.CalculatedAirdensity_GRe(IsConsidered));
TurbineData(ii-2).HWShub1 = mean(newdata.HWShub1(IsConsidered));
TurbineData(ii-2).av100msabswinddirectionGRe = mean(newdata.av100msabswinddirectionGRe(IsConsidered));
end
The function ReadFamosDataIntoTimeTable is to convert the .dat binary data into .mat dataI am working on a project where I have to analyse high frequency wind data (50 Hz) coming from a wind turbine data measurement system. this amount of data must be converted from .dat binary files to .mat files which I can use in matlab. The data must then be filtered and then averaged over 10 minutes to be compared to the data from another measurement system. Doing all this requires the analysis of thousands of data and it’s very time consuming (right now it is about 105 seconds just for the data of 2 days). How can I speed the process up?
for ii = 3:length(filedir)
filename = filedir(ii).name;
newdata = ReadFamosDataIntoTimeTable(filename);
% filter
IsConsidered = newdata.avbladeangleGRe<40 … % normal operation
& newdata.RAWS9>0.5 … % good availability
& ~isnan(newdata.iv10mswindspeed2GRe) & newdata.av100msabswinddirectionGRe>180 & newdata.av100msabswinddirectionGRe<250;
TurbineData(ii-2).date = filename;
TurbineData(ii-2).iv10mswindspeed2GRe = mean(newdata.iv10mswindspeed2GRe(IsConsidered));
TurbineData(ii-2).ivactivepowerGRe = mean(newdata.ivactivepowerGRe(IsConsidered));
TurbineData(ii-2).CalculatedAirdensity_GRe = mean(newdata.CalculatedAirdensity_GRe(IsConsidered));
TurbineData(ii-2).HWShub1 = mean(newdata.HWShub1(IsConsidered));
TurbineData(ii-2).av100msabswinddirectionGRe = mean(newdata.av100msabswinddirectionGRe(IsConsidered));
end
The function ReadFamosDataIntoTimeTable is to convert the .dat binary data into .mat data I am working on a project where I have to analyse high frequency wind data (50 Hz) coming from a wind turbine data measurement system. this amount of data must be converted from .dat binary files to .mat files which I can use in matlab. The data must then be filtered and then averaged over 10 minutes to be compared to the data from another measurement system. Doing all this requires the analysis of thousands of data and it’s very time consuming (right now it is about 105 seconds just for the data of 2 days). How can I speed the process up?
for ii = 3:length(filedir)
filename = filedir(ii).name;
newdata = ReadFamosDataIntoTimeTable(filename);
% filter
IsConsidered = newdata.avbladeangleGRe<40 … % normal operation
& newdata.RAWS9>0.5 … % good availability
& ~isnan(newdata.iv10mswindspeed2GRe) & newdata.av100msabswinddirectionGRe>180 & newdata.av100msabswinddirectionGRe<250;
TurbineData(ii-2).date = filename;
TurbineData(ii-2).iv10mswindspeed2GRe = mean(newdata.iv10mswindspeed2GRe(IsConsidered));
TurbineData(ii-2).ivactivepowerGRe = mean(newdata.ivactivepowerGRe(IsConsidered));
TurbineData(ii-2).CalculatedAirdensity_GRe = mean(newdata.CalculatedAirdensity_GRe(IsConsidered));
TurbineData(ii-2).HWShub1 = mean(newdata.HWShub1(IsConsidered));
TurbineData(ii-2).av100msabswinddirectionGRe = mean(newdata.av100msabswinddirectionGRe(IsConsidered));
end
The function ReadFamosDataIntoTimeTable is to convert the .dat binary data into .mat data optimization, for loop, speed, data conversion MATLAB Answers — New Questions

​

Hi,
Say I have a 5 by 2 matrix in the form:
A = [2 5; 9 7; 10 2; 3 2; 1 9]
And I want to make a 10 by 1 matrix from these values so that I get:
B = [2;5;9;7;10;2;3;2;1;9]
How would I do this?

I know there is a probably a simple fix, but I haven’t been able to do it.
Many thanks,

ScottHi,
Say I have a 5 by 2 matrix in the form:
A = [2 5; 9 7; 10 2; 3 2; 1 9]
And I want to make a 10 by 1 matrix from these values so that I get:
B = [2;5;9;7;10;2;3;2;1;9]
How would I do this?

I know there is a probably a simple fix, but I haven’t been able to do it.
Many thanks,

Scott Hi,
Say I have a 5 by 2 matrix in the form:
A = [2 5; 9 7; 10 2; 3 2; 1 9]
And I want to make a 10 by 1 matrix from these values so that I get:
B = [2;5;9;7;10;2;3;2;1;9]
How would I do this?

I know there is a probably a simple fix, but I haven’t been able to do it.
Many thanks,

Scott matrices, matrix manipulation MATLAB Answers — New Questions

​

I’m working on a numerical model involving the merging of four turbulent plumes, each originating from a circular area source located at the corners of a square. My goal is to compute and visualize the merged outer contour that encloses all four interacting plumes. I have attached the code below, the equation I am solving is at the bottom of this code. When I make the plot some gaps appear and is there a better way to plot it so the gaps disappear
%% Four-Plume Contour + Area Plotting with Automatic k_crit Detection
close all; clc; clf
global r0 k theta
%%Parameter Input
r0 = 0.25; %Plume Source Radius
%% STEP 1: Find k_crit by setting θ = π/4 and solving corresponding function, g(rho) = 0 for r0
theta = pi/4;
kL = 0*(1 – r0^2)^2; % Conservative lower bound
kU = 1.05*(1 – r0^2)^2; % Previous Theoretical maximum
k_testArray = linspace(kL,kU, 501);
rho_valid = NaN(size(k_testArray)); % Store valid roots for each k in sweep
k_crit = NaN; % First k where g(rho)=0 has a real root
contactDetected = false; % Logical flag to stop at first contact
for jj = 1:length(k_testArray)
k = k_testArray(jj);
rho = mnhf_secant(@FourPlumes_Equation, [0.1 0.5], 1e-6, 0);
if rho > 0 && isreal(rho) && isfinite(rho)
rho_valid(jj) = rho;
if ~contactDetected
k_crit = k;
contactDetected = true;
end
end
end
fprintf(‘Detected critical k value (first contact with θ = π/4): k_crit = %.8fnn’, k_crit);
%% Parameter input: k values
Nt = 10001; % Resolution, must be an odd number
theta_array = linspace(-pi/2, pi/2, Nt);
k_array = [0.6]
%k_array = [0.2 0.3 0.4 0.5 0.6 0.7 k_crit 0.8 0.9 1.2 1.5]; % test range values for k
area = zeros(1,length(k_array));
%% STEP 2: Loop over k values, plot, and compute area
figure(1); hold on; box on
for jj = 1:length(k_array)
k = k_array(jj);
% For merged Case
if k >= k_crit
rho_array1 = zeros(1, Nt);
for ii = 1:Nt
theta = theta_array(ii);
rho_array1(ii) = mnhf_secant(@FourPlumes_Equation, [1.5 2], 1e-6, 0);
if rho_array1(ii) < 0
rho_array1(ii) = NaN;
end
end
[x1, y1] = pol2cart(theta_array, rho_array1);
plot(x1, y1, ‘-k’); axis square; xlim([0 2]); ylim([-1 1]);
% Plot θ = π/4 lines for visual reference
x_diag = linspace(0, 2, Nt); y_diag = x_diag;
plot(x_diag, y_diag, ‘r–‘); plot(x_diag, -y_diag, ‘r–‘);
% Area estimation using trapezoidal rule
iq = find(y1 >= y_diag); % intersection
if ~isempty(iq)
area(jj) = 2 * abs(trapz(x1(min(iq):(Nt-1)/2+1), y1(min(iq):(Nt-1)/2+1)));
area(jj) = area(jj) + x1(min(iq)-1)^2;
end
else
% Before Merging Case
rho_array1 = zeros(1, Nt);
for ii = 1:Nt
theta = theta_array(ii);
rho_array1(ii) = mnhf_secant(@FourPlumes_Equation, [1.5 2], 1e-6, 0);
if rho_array1(ii) < 0
rho_array1(ii) = NaN;
end
end
[x1, y1] = pol2cart(theta_array, rho_array1);
plot(x1, y1, ‘-k’); axis square;
% Detect breakdown region and recompute in "gap"
ij = find(isnan(rho_array1(1:round(end/2))));
if ~isempty(ij)
theta_min = theta_array(max(ij));
theta_array2 = linspace(theta_min, -theta_min, Nt);
rho_array2 = zeros(1, Nt);
for ii = 1:Nt
theta = theta_array2(ii);
rho_array2(ii) = mnhf_secant(@FourPlumes_Equation, [0.1 0.9], 1e-6, 0);
if rho_array2(ii) < 0
rho_array2(ii) = NaN;
end
end
[x2, y2] = pol2cart(theta_array2, rho_array2);
plot(x2, y2, ‘r.’);
ik = (Nt-1)/2 + find(isnan(x2((Nt-1)/2+1:end)));
x = [x2((Nt-1)/2+1:min(ik)-1), x1(max(ij)+1:(Nt-1)/2+1)];
y = [y2((Nt-1)/2+1:min(ik)-1), -y1(max(ij)+1:(Nt-1)/2+1)];
plot(x, y, ‘b–‘); plot(x, -y, ‘b–‘);
% Area from symmetric region
area(jj) = 2 * trapz(x, y);
end
end
% Plot plume source circles (n = 4)
pos1 = [1 – r0, -r0, 2*r0, 2*r0]; % (1, 0)
pos2 = [-1 – r0, -r0, 2*r0, 2*r0]; % (-1, 0)
pos3 = [-r0, 1 – r0, 2*r0, 2*r0]; % (0, 1)
pos4 = [-r0, -1 – r0, 2*r0, 2*r0]; % (0, -1)
rectangle(‘Position’, pos1, ‘Curvature’, [1 1], ‘FaceColor’, ‘k’);
rectangle(‘Position’, pos2, ‘Curvature’, [1 1], ‘FaceColor’, ‘k’);
rectangle(‘Position’, pos3, ‘Curvature’, [1 1], ‘FaceColor’, ‘k’);
rectangle(‘Position’, pos4, ‘Curvature’, [1 1], ‘FaceColor’, ‘k’);
end
%% Function: Four-Plume Equation (Li & Flynn B3) —
function f = FourPlumes_Equation(rho)
global r0 k theta
f = (rho.^2 – 2*rho*cos(theta) + 1 – r0^2) .* …
(rho.^2 – 2*rho*cos(theta + pi/2) + 1 – r0^2) .* …
(rho.^2 – 2*rho*cos(theta + pi) + 1 – r0^2) .* …
(rho.^2 – 2*rho*cos(theta + 1.5*pi) + 1 – r0^2) – k^2;
end

function r=mnhf_secant(Fun,x,tol,trace)
%MNHF_SECANT finds the root of "Fun" using secant scheme.
%
% Fun – name of the external function
% x – vector of length 2, (initial guesses)
% tol – tolerance
% trace – print intermediate results
%
% Usage mnhf_secant(@poly1,[-0.5 2.0],1e-8,1)
% poly1 is the name of the external function.
% [-0.5 2.0] are the initial guesses for the root.

% Check inputs.
if nargin < 4, trace = 1; end
if nargin < 3, tol = 1e-8; end
if (length(x) ~= 2)
error(‘Please provide two initial guesses’)
end

f = feval(Fun,x); % Fun is assumed to accept a vector

for i = 1:100
x3 = x(1)-f(1)*(x(2)-x(1))/(f(2)-f(1)); % Update the guess.
f3 = feval(Fun,x3); % Function evaluation.
if trace, fprintf(1,’%3i %12.5f %12.5fn’, i,x3,f3); end
if abs(f3) < tol % Check for convergenece.
r = x3;
return
else % Reset values for x(1), x(2), f(1) and f(2).
x(1) = x(2); f(1) = f(2); x(2) = x3; f(2) = f3;
end
end

r = NaN;
endI’m working on a numerical model involving the merging of four turbulent plumes, each originating from a circular area source located at the corners of a square. My goal is to compute and visualize the merged outer contour that encloses all four interacting plumes. I have attached the code below, the equation I am solving is at the bottom of this code. When I make the plot some gaps appear and is there a better way to plot it so the gaps disappear
%% Four-Plume Contour + Area Plotting with Automatic k_crit Detection
close all; clc; clf
global r0 k theta
%%Parameter Input
r0 = 0.25; %Plume Source Radius
%% STEP 1: Find k_crit by setting θ = π/4 and solving corresponding function, g(rho) = 0 for r0
theta = pi/4;
kL = 0*(1 – r0^2)^2; % Conservative lower bound
kU = 1.05*(1 – r0^2)^2; % Previous Theoretical maximum
k_testArray = linspace(kL,kU, 501);
rho_valid = NaN(size(k_testArray)); % Store valid roots for each k in sweep
k_crit = NaN; % First k where g(rho)=0 has a real root
contactDetected = false; % Logical flag to stop at first contact
for jj = 1:length(k_testArray)
k = k_testArray(jj);
rho = mnhf_secant(@FourPlumes_Equation, [0.1 0.5], 1e-6, 0);
if rho > 0 && isreal(rho) && isfinite(rho)
rho_valid(jj) = rho;
if ~contactDetected
k_crit = k;
contactDetected = true;
end
end
end
fprintf(‘Detected critical k value (first contact with θ = π/4): k_crit = %.8fnn’, k_crit);
%% Parameter input: k values
Nt = 10001; % Resolution, must be an odd number
theta_array = linspace(-pi/2, pi/2, Nt);
k_array = [0.6]
%k_array = [0.2 0.3 0.4 0.5 0.6 0.7 k_crit 0.8 0.9 1.2 1.5]; % test range values for k
area = zeros(1,length(k_array));
%% STEP 2: Loop over k values, plot, and compute area
figure(1); hold on; box on
for jj = 1:length(k_array)
k = k_array(jj);
% For merged Case
if k >= k_crit
rho_array1 = zeros(1, Nt);
for ii = 1:Nt
theta = theta_array(ii);
rho_array1(ii) = mnhf_secant(@FourPlumes_Equation, [1.5 2], 1e-6, 0);
if rho_array1(ii) < 0
rho_array1(ii) = NaN;
end
end
[x1, y1] = pol2cart(theta_array, rho_array1);
plot(x1, y1, ‘-k’); axis square; xlim([0 2]); ylim([-1 1]);
% Plot θ = π/4 lines for visual reference
x_diag = linspace(0, 2, Nt); y_diag = x_diag;
plot(x_diag, y_diag, ‘r–‘); plot(x_diag, -y_diag, ‘r–‘);
% Area estimation using trapezoidal rule
iq = find(y1 >= y_diag); % intersection
if ~isempty(iq)
area(jj) = 2 * abs(trapz(x1(min(iq):(Nt-1)/2+1), y1(min(iq):(Nt-1)/2+1)));
area(jj) = area(jj) + x1(min(iq)-1)^2;
end
else
% Before Merging Case
rho_array1 = zeros(1, Nt);
for ii = 1:Nt
theta = theta_array(ii);
rho_array1(ii) = mnhf_secant(@FourPlumes_Equation, [1.5 2], 1e-6, 0);
if rho_array1(ii) < 0
rho_array1(ii) = NaN;
end
end
[x1, y1] = pol2cart(theta_array, rho_array1);
plot(x1, y1, ‘-k’); axis square;
% Detect breakdown region and recompute in "gap"
ij = find(isnan(rho_array1(1:round(end/2))));
if ~isempty(ij)
theta_min = theta_array(max(ij));
theta_array2 = linspace(theta_min, -theta_min, Nt);
rho_array2 = zeros(1, Nt);
for ii = 1:Nt
theta = theta_array2(ii);
rho_array2(ii) = mnhf_secant(@FourPlumes_Equation, [0.1 0.9], 1e-6, 0);
if rho_array2(ii) < 0
rho_array2(ii) = NaN;
end
end
[x2, y2] = pol2cart(theta_array2, rho_array2);
plot(x2, y2, ‘r.’);
ik = (Nt-1)/2 + find(isnan(x2((Nt-1)/2+1:end)));
x = [x2((Nt-1)/2+1:min(ik)-1), x1(max(ij)+1:(Nt-1)/2+1)];
y = [y2((Nt-1)/2+1:min(ik)-1), -y1(max(ij)+1:(Nt-1)/2+1)];
plot(x, y, ‘b–‘); plot(x, -y, ‘b–‘);
% Area from symmetric region
area(jj) = 2 * trapz(x, y);
end
end
% Plot plume source circles (n = 4)
pos1 = [1 – r0, -r0, 2*r0, 2*r0]; % (1, 0)
pos2 = [-1 – r0, -r0, 2*r0, 2*r0]; % (-1, 0)
pos3 = [-r0, 1 – r0, 2*r0, 2*r0]; % (0, 1)
pos4 = [-r0, -1 – r0, 2*r0, 2*r0]; % (0, -1)
rectangle(‘Position’, pos1, ‘Curvature’, [1 1], ‘FaceColor’, ‘k’);
rectangle(‘Position’, pos2, ‘Curvature’, [1 1], ‘FaceColor’, ‘k’);
rectangle(‘Position’, pos3, ‘Curvature’, [1 1], ‘FaceColor’, ‘k’);
rectangle(‘Position’, pos4, ‘Curvature’, [1 1], ‘FaceColor’, ‘k’);
end
%% Function: Four-Plume Equation (Li & Flynn B3) —
function f = FourPlumes_Equation(rho)
global r0 k theta
f = (rho.^2 – 2*rho*cos(theta) + 1 – r0^2) .* …
(rho.^2 – 2*rho*cos(theta + pi/2) + 1 – r0^2) .* …
(rho.^2 – 2*rho*cos(theta + pi) + 1 – r0^2) .* …
(rho.^2 – 2*rho*cos(theta + 1.5*pi) + 1 – r0^2) – k^2;
end

function r=mnhf_secant(Fun,x,tol,trace)
%MNHF_SECANT finds the root of "Fun" using secant scheme.
%
% Fun – name of the external function
% x – vector of length 2, (initial guesses)
% tol – tolerance
% trace – print intermediate results
%
% Usage mnhf_secant(@poly1,[-0.5 2.0],1e-8,1)
% poly1 is the name of the external function.
% [-0.5 2.0] are the initial guesses for the root.

% Check inputs.
if nargin < 4, trace = 1; end
if nargin < 3, tol = 1e-8; end
if (length(x) ~= 2)
error(‘Please provide two initial guesses’)
end

f = feval(Fun,x); % Fun is assumed to accept a vector

for i = 1:100
x3 = x(1)-f(1)*(x(2)-x(1))/(f(2)-f(1)); % Update the guess.
f3 = feval(Fun,x3); % Function evaluation.
if trace, fprintf(1,’%3i %12.5f %12.5fn’, i,x3,f3); end
if abs(f3) < tol % Check for convergenece.
r = x3;
return
else % Reset values for x(1), x(2), f(1) and f(2).
x(1) = x(2); f(1) = f(2); x(2) = x3; f(2) = f3;
end
end

r = NaN;
end I’m working on a numerical model involving the merging of four turbulent plumes, each originating from a circular area source located at the corners of a square. My goal is to compute and visualize the merged outer contour that encloses all four interacting plumes. I have attached the code below, the equation I am solving is at the bottom of this code. When I make the plot some gaps appear and is there a better way to plot it so the gaps disappear
%% Four-Plume Contour + Area Plotting with Automatic k_crit Detection
close all; clc; clf
global r0 k theta
%%Parameter Input
r0 = 0.25; %Plume Source Radius
%% STEP 1: Find k_crit by setting θ = π/4 and solving corresponding function, g(rho) = 0 for r0
theta = pi/4;
kL = 0*(1 – r0^2)^2; % Conservative lower bound
kU = 1.05*(1 – r0^2)^2; % Previous Theoretical maximum
k_testArray = linspace(kL,kU, 501);
rho_valid = NaN(size(k_testArray)); % Store valid roots for each k in sweep
k_crit = NaN; % First k where g(rho)=0 has a real root
contactDetected = false; % Logical flag to stop at first contact
for jj = 1:length(k_testArray)
k = k_testArray(jj);
rho = mnhf_secant(@FourPlumes_Equation, [0.1 0.5], 1e-6, 0);
if rho > 0 && isreal(rho) && isfinite(rho)
rho_valid(jj) = rho;
if ~contactDetected
k_crit = k;
contactDetected = true;
end
end
end
fprintf(‘Detected critical k value (first contact with θ = π/4): k_crit = %.8fnn’, k_crit);
%% Parameter input: k values
Nt = 10001; % Resolution, must be an odd number
theta_array = linspace(-pi/2, pi/2, Nt);
k_array = [0.6]
%k_array = [0.2 0.3 0.4 0.5 0.6 0.7 k_crit 0.8 0.9 1.2 1.5]; % test range values for k
area = zeros(1,length(k_array));
%% STEP 2: Loop over k values, plot, and compute area
figure(1); hold on; box on
for jj = 1:length(k_array)
k = k_array(jj);
% For merged Case
if k >= k_crit
rho_array1 = zeros(1, Nt);
for ii = 1:Nt
theta = theta_array(ii);
rho_array1(ii) = mnhf_secant(@FourPlumes_Equation, [1.5 2], 1e-6, 0);
if rho_array1(ii) < 0
rho_array1(ii) = NaN;
end
end
[x1, y1] = pol2cart(theta_array, rho_array1);
plot(x1, y1, ‘-k’); axis square; xlim([0 2]); ylim([-1 1]);
% Plot θ = π/4 lines for visual reference
x_diag = linspace(0, 2, Nt); y_diag = x_diag;
plot(x_diag, y_diag, ‘r–‘); plot(x_diag, -y_diag, ‘r–‘);
% Area estimation using trapezoidal rule
iq = find(y1 >= y_diag); % intersection
if ~isempty(iq)
area(jj) = 2 * abs(trapz(x1(min(iq):(Nt-1)/2+1), y1(min(iq):(Nt-1)/2+1)));
area(jj) = area(jj) + x1(min(iq)-1)^2;
end
else
% Before Merging Case
rho_array1 = zeros(1, Nt);
for ii = 1:Nt
theta = theta_array(ii);
rho_array1(ii) = mnhf_secant(@FourPlumes_Equation, [1.5 2], 1e-6, 0);
if rho_array1(ii) < 0
rho_array1(ii) = NaN;
end
end
[x1, y1] = pol2cart(theta_array, rho_array1);
plot(x1, y1, ‘-k’); axis square;
% Detect breakdown region and recompute in "gap"
ij = find(isnan(rho_array1(1:round(end/2))));
if ~isempty(ij)
theta_min = theta_array(max(ij));
theta_array2 = linspace(theta_min, -theta_min, Nt);
rho_array2 = zeros(1, Nt);
for ii = 1:Nt
theta = theta_array2(ii);
rho_array2(ii) = mnhf_secant(@FourPlumes_Equation, [0.1 0.9], 1e-6, 0);
if rho_array2(ii) < 0
rho_array2(ii) = NaN;
end
end
[x2, y2] = pol2cart(theta_array2, rho_array2);
plot(x2, y2, ‘r.’);
ik = (Nt-1)/2 + find(isnan(x2((Nt-1)/2+1:end)));
x = [x2((Nt-1)/2+1:min(ik)-1), x1(max(ij)+1:(Nt-1)/2+1)];
y = [y2((Nt-1)/2+1:min(ik)-1), -y1(max(ij)+1:(Nt-1)/2+1)];
plot(x, y, ‘b–‘); plot(x, -y, ‘b–‘);
% Area from symmetric region
area(jj) = 2 * trapz(x, y);
end
end
% Plot plume source circles (n = 4)
pos1 = [1 – r0, -r0, 2*r0, 2*r0]; % (1, 0)
pos2 = [-1 – r0, -r0, 2*r0, 2*r0]; % (-1, 0)
pos3 = [-r0, 1 – r0, 2*r0, 2*r0]; % (0, 1)
pos4 = [-r0, -1 – r0, 2*r0, 2*r0]; % (0, -1)
rectangle(‘Position’, pos1, ‘Curvature’, [1 1], ‘FaceColor’, ‘k’);
rectangle(‘Position’, pos2, ‘Curvature’, [1 1], ‘FaceColor’, ‘k’);
rectangle(‘Position’, pos3, ‘Curvature’, [1 1], ‘FaceColor’, ‘k’);
rectangle(‘Position’, pos4, ‘Curvature’, [1 1], ‘FaceColor’, ‘k’);
end
%% Function: Four-Plume Equation (Li & Flynn B3) —
function f = FourPlumes_Equation(rho)
global r0 k theta
f = (rho.^2 – 2*rho*cos(theta) + 1 – r0^2) .* …
(rho.^2 – 2*rho*cos(theta + pi/2) + 1 – r0^2) .* …
(rho.^2 – 2*rho*cos(theta + pi) + 1 – r0^2) .* …
(rho.^2 – 2*rho*cos(theta + 1.5*pi) + 1 – r0^2) – k^2;
end

function r=mnhf_secant(Fun,x,tol,trace)
%MNHF_SECANT finds the root of "Fun" using secant scheme.
%
% Fun – name of the external function
% x – vector of length 2, (initial guesses)
% tol – tolerance
% trace – print intermediate results
%
% Usage mnhf_secant(@poly1,[-0.5 2.0],1e-8,1)
% poly1 is the name of the external function.
% [-0.5 2.0] are the initial guesses for the root.

% Check inputs.
if nargin < 4, trace = 1; end
if nargin < 3, tol = 1e-8; end
if (length(x) ~= 2)
error(‘Please provide two initial guesses’)
end

f = feval(Fun,x); % Fun is assumed to accept a vector

for i = 1:100
x3 = x(1)-f(1)*(x(2)-x(1))/(f(2)-f(1)); % Update the guess.
f3 = feval(Fun,x3); % Function evaluation.
if trace, fprintf(1,’%3i %12.5f %12.5fn’, i,x3,f3); end
if abs(f3) < tol % Check for convergenece.
r = x3;
return
else % Reset values for x(1), x(2), f(1) and f(2).
x(1) = x(2); f(1) = f(2); x(2) = x3; f(2) = f3;
end
end

r = NaN;
end numerical solution, nonlinear-equations, root-finding, polar-coordinates, contour-gaps MATLAB Answers — New Questions

​