numerical Instabilities for bessel functions

306 vues (au cours des 30 derniers jours)
Javeria
Javeria le 8 Déc 2025
Commenté : Torsten le 29 Jan 2026 à 13:50
I write the following code but i used the following parametrs and i did not get any numerical instabilities:
RI = 34.5; % Inner radius (m)
ratio = 47.5/34.5; % Ratio R_E / R_I
RE = ratio*RI; % outer radius (m)
h = 200/34.5 * RI; % Water depth (m)
d = 38/34.5 * RI; % Interface depth (m)
b = h-d;
But when i switch the parameters with this then i get numerical instabilities
RI = 0.4; % Inner radius (m)
RE = 0.8; % Outer radius (m)
d = 0.3; % Draft (m)
h = 1.0; % Water depth (m)
b = h-d; %from cylinder bottom to seabed (m)
the following is code
clear all;
close all;
clc;
tic;
% diary('iter_log.txt');
% diary on
% fprintf('Run started: %s\n', datestr(now));
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Parameters from your setup
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
N = 10; % Number of N modes
M = 10; % Number of M modes
Q = 10; % Number of Q modes
RI = 0.4; % Inner radius (m)
RE = 0.8; % Outer radius (m)
d = 0.3; % Draft (m)
h = 1.0; % Water depth (m) % Interface depth (m)
b = h-d; % Distance from cylinder bottom to seabed (m)
g = 9.81; % Gravity (m/s^2)
tau =0.2; % porosity ratio
gamma = 1.0; % discharge coefficient
b1 = (1 - tau) / (2 * gamma * tau^2); % nonlinear coeff
tol = 1e-4; % convergence tolerance
max_iter = 100; % max iterations
l = 0; % Azimuthal order
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% === Compute Bessel roots and scaled chi values ===
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
roots = bessel0j(l, Q); % Zeros of J_l(x)
chi = roots ./ RI; % chi_q^l = roots scaled by radius
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% ===== Define Range for k_0
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% X_kRE = linspace(0.0025, 0.16,100); %%%%%% this is now k_0
X_kRE = linspace(0.05, 4.5, 40);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% ==============Initialize eta0 array before loop
%===============Initialize omega array before loop
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
eta0 = zeros(size(X_kRE)); % Preallocate
omega = zeros(size(X_kRE)); %%%% omega preallocate
% iters_used = zeros(length(X_kRE),1);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for idx = 1:length(X_kRE)
wk = X_kRE(idx); % dimensional wavenumber
omega(idx) = sqrt(g * wk * tanh(wk * h)); % dispersion relation
% fprintf('\n--- idx %3d (T=%6.3f s) ---\n', idx, 2*pi/omega(idx));
% drawnow;
%%%%%%%%%%%%% Compute group velocity %%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Cg(idx) = (g * tanh(wk * h) + g * wk * h * (sech(wk * h))^2) ...
* omega(idx) / (2 * g * wk * tanh(wk * h));
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% =======Compute a1 based on tau
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
a1(idx) = 0.93 * (1 - tau) * Cg(idx);
% dissip(idx) = a1(idx)/omega(idx);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%=============== Find derivative of Z0 at z = -d
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Z0d = (cosh(wk * h) * wk * sinh(wk * b)) / (2 * wk * h + sinh(2 * wk * h));
fun_alpha = @(x) omega(idx)^2 + g*x.*tan(x*h); %dispersion equation for evanescent modes
opts_lsq = optimoptions('lsqnonlin','Display','off'); % ← add this
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,opts_lsq);
%%%%%%%%%%%%% Derivative of Z_n at z = -d %%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Znd(n) = -k(n) * (cos(k(n)*h) * sin(k(n)*b)) / (2*k(n)*h + sin(2*k(n)*h));
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
% % %%%%%%%%%% Write Matrix W %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
W = zeros(N,Q);
% %%W_{0q}
for q = 1:Q
r = exp(-2*chi(q)*b); % e^{-2*chi*b}
coth_chib = (1 + r) / (1 - r); % coth(chi*b) via decaying exp
W(1,q) = -(4*chi(q)/cosh(wk*h))*(( wk*sinh(wk*b)*coth_chib - chi(q)*cosh(wk*b) ) / ( wk^2 - chi(q)^2 ) ); %prof.chen
%
% % W(1, q) = - (4 * chi(q) / cosh(wk * h)) * (chi(q) * sinh(chi(q) * b) * cosh(wk * b) - wk * sinh(wk * b) * cosh(chi(q) * b)) ...
% % / ((chi(q)^2 - wk^2) * sinh(chi(q) * b));
end
% % %%% Write W_{nq}
for n = 2:N
for q = 1:Q
r = exp(-2*chi(q)*b); % e^{-2*chi*b}
coth_chib = (1 + r) / (1 - r); % coth(chi*b), stable
W(n,q) = -(4*chi(q)/cos(k(n)*h)) * ...
( ( k(n)*sin(k(n)*b)*coth_chib + chi(q)*cos(k(n)*b) ) ... %prof.chen
/ ( k(n)^2 + chi(q)^2 ) );
% % W(n, q) = - (4 * chi(q) / cos(k(n) * h)) * (chi(q) * sinh(chi(q) * b) * cos(k(n) * b) + k(n) * sin(k(n) * b) * cosh(chi(q) * b)) ...
% % / ((chi(q)^2 + k(n)^2) * sinh(chi(q) * b));
end
end
%%%%%%%%%%%%%%%%%%%%%%%Find E1_{q}%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
f2 = omega(idx)^2 / g;
E1 = zeros(Q,1);
S_q = zeros(Q,1);
for q = 1:Q
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% this expression is same as prof. Chen code
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
exp2chi_d = exp(-2 * chi(q) * d); % e^{-2*chi*d}
exp2chi_b = exp(-2 * chi(q) * b); % e^{-2*chi*(h-d)}
exp2chi_h = exp(-2 * chi(q) * h); % e^{-2*chi*h}
exp1chi_d = exp(-chi(q) * d); % e^{-chi*d}
% intermediate ratio (same order as Fortran)
BB = (f2 - chi(q)) * (exp2chi_b / exp2chi_d);
BB = BB + (f2 + chi(q)) * (1 + 2 * exp2chi_b);
BB = BB / (f2 - chi(q) + (f2 + chi(q)) * exp2chi_d);
BB = BB / (1 + exp2chi_b);
% final coefficient
E1(q) = 2 * (BB * exp2chi_d - 1);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% This expression is that one derived directly
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% % N0, D0, H0 exactly as in the math above
% N0 = (chi(q) - f2) + (chi(q) + f2) * exp2chi_h;
% D0 = (f2 - chi(q)) + (f2 + chi(q)) * exp2chi_d;
% H0 = 1 + exp2chi_b;
%
% % final coefficient (no add/subtract step)
% E1(q) = 2 * N0 / (D0 * H0);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% This one is the hyperbolic one form
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% num = chi(q)*cosh(chi(q)*h) - f2*sinh(chi(q)*h);
% den = (f2*cosh(chi(q)*d) - chi(q)*sinh(chi(q)*d)) * cosh(chi(q)*(h - d)); % or: * cosh(chi(q)*b)
% E1(q) = num / den;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%% This is the upper region velocity potential at z=0, write
%%%%%%%%%%%%%%%% from prof.chen code
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
CCmc = 1 + (exp2chi_b / exp2chi_d) + 2*exp2chi_b;
CCmc = CCmc / (1 + exp2chi_b);
CCmc = BB * (1 + exp2chi_d) - CCmc;
S_q(q) = CCmc * exp1chi_d; % this equals the summand value for this q
% --- direct hyperbolic (exact) ---
% % % C_hyp = ( chi(q)*cosh(chi(q)*h) - f2*sinh(chi(q)*h) ) ...
% % % / ( (f2*cosh(chi(q)*d) - chi(q)*sinh(chi(q)*d)) * cosh(chi(q)*b) );
% % %
% % % S_hyp(q) = C_hyp * cosh(chi(q)*d) + sinh(chi(q)*h) / cosh(chi(q)*b); % summand
end
% H(q): diagonal term (exponential form)
H1 = zeros(Q,1);
for q = 1:Q
% dissip =
r = exp(-2*chi(q)*b); % e^{-2*chi*b}
mid = 4*r/(1 - r^2) - E1(q); % 2/sinh(2*chi*b) - C_m
H1(q) = mid + 1i*a1(idx)*chi(q)/omega(idx); % mid + i*(a1/w)*chi
end
H = diag(H1);
%%%%%%%%%%matric G%%%%%%%%%%%%%%%%%%
G = zeros(Q, N); % Preallocate
for q = 1:Q
G(q,1) = 1i * 2*(a1(idx)/(omega(idx)*RI))* Z0d ...
* ( chi(q) / (chi(q)^2 - wk^2) )* ( besselj(l, wk*RI) / dbesselj(l, wk*RI) ); % this one is my calculated
end
% %G_qn
for q = 1:Q
for n = 2:N
G(q, n) = 1i * 2*(a1(idx)/(omega(idx)*RI))* Znd(n) ...
* ( chi(q) / (k(n)^2 + chi(q)^2) )* ( besseli(l, k(n)*RI) / dbesseli(l, k(n)*RI) ); % this one is my calculated
end
end
% % % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% %Define the right hand side vector
U = zeros(2*M + 2*N + Q, 1); % Full vector, M+N+M+N+Q elements
% % Block 1: U (size M = 4)
for i = 1:M
if i == 1
U(i, 1) = (sinh(wk*b))/((b*wk)*cosh(wk*h))*besselj(l, wk*RE) ; % Z_0^l
else
U(i, 1) = (2 *wk*(-1)^(i-1)*sinh(wk*b))/(b *(wk^2 +m(i)^2)*cosh(wk*h))*besselj(l, wk*RE); %Z_m^l
end
end
% % Block 2: Y (size N = 3)
for j = 1:N
if j == 1
U(j + M, 1) = -dbesselj(l, wk*RE) * (2*wk*h + sinh(2*wk*h)) /(cosh(wk*h)^2); % Y_0^l
else
U(j + M, 1) = 0; % Y_n^l
end
end
% % Block 3: X (size M = 4)
for i = 1:M
U(i + M + N, 1) = 0; % X_0^l, X_m^l
end
% Block 4: W (size N = 3)
for j = 1:N
U(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);
ZMQ = zeros(M,Q);
ZNQ = zeros(N,Q);
ZQM = zeros(Q,M);
ZQN = zeros(Q,N);
% Construct the full block matrix S
S = [ I_M, -A, ZMM, ZMN, ZMQ;
-B, I_N, -C, ZNN, ZNQ;
ZMM, ZMN, I_M, -D, ZMQ;
-E, ZNN, -F, I_N, -W;
ZQM, ZQN, ZQM, -G, H ];
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
converged = false;
psi = zeros(Q,1);
T_old = []; % <-- track the full solution from prev. iter
% T_old = zeros(2*M + 2*N + Q, 1); % previous unknowns (seed)
for iter = 1:max_iter
% (A) overwrite ONLY Block 5 of U from current psi
for q = 1:Q
U(2*M + 2*N + q) = -(1i*b1/omega(idx)) * (2/(RI^2 * dbesselj(l, chi(q)*RI))) * psi(q);
end
T = S \ U; % Solve the linear system
% T = pinv(S) * U; % Use pseudoinverse for stability
b_vec = T(1:M); % Coefficients b^l
a_vec = T(M+1 : M+N); % Coefficients a^l
c_vec = T(M+N+1 : 2*M+N); % Coefficients c^l
d_vec = T(2*M+N+1:2*M+2*N); % Coefficients d^l
e_vec = T(2*M+2*N+1:end); % (Q×1)
% (D) update psi for NEXT iteration from CURRENT coefficients
for q = 1:Q
integrand = @(r) abs(v_D(N,Q,r,Z0d,wk,RI,l,d_vec,Znd,k,e_vec,chi)) ...
.* v_D(N,Q,r,Z0d,wk,RI,l,d_vec,Znd,k,e_vec,chi) ...
.* besselj(l, chi(q)*r) .* r;
psi(q) = integral(integrand, 0, RI, 'AbsTol',1e-8, 'RelTol',1e-6);
end
% % % drawnow; % optional: flush output each iter
% === compact per-iteration print, only for T in [5,35] ===
% Tcur = 2*pi/omega(idx);
% if Tcur >= 5 && Tcur <= 35
% fprintf('idx %3d | iter %2d | ||psi|| = %.3e\n', idx, iter, norm(psi));
% drawnow;
% end
% % % % (E) convergence on ALL unknowns (T vs previous T)
% --- per-iteration diagnostics (optional) ---
if ~isempty(T_old)
diff_T = max(abs(T - T_old));
% Tcur = 2*pi/omega(idx);
% if Tcur >= 5 && Tcur <= 35
% fprintf('k0 idx %3d | iter %2d: max|ΔT| = %.3e, ||psi|| = %.3e\n', ...
% idx, iter, diff_T, norm(psi));
% drawnow; % ← add this to flush each iteration line too
% end
if diff_T < tol
converged = true;
break
end
end
T_old = T;
end % <<< end of iter loop
% % ---- ONE summary line per frequency (place it here) ----
% iters_used(idx) = iter; % how many iterations this frequency needed
% fprintf('idx %3d (T=%6.3f s): iters=%2d\n', idx, 2*pi/omega(idx), iter);
% drawnow;
% -------------- end nonlinear iteration --------------
% % %%%%%%%%% 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
phiUell = 0;
for q =1:Q
phiUell = phiUell +e_vec(q)* S_q(q);
end
eta0(idx) = abs(term1+sum1+phiUell);
end
% % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% % % Plotting Data
% % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% %
% %
% figure(2); clf;
% >>> CHANGES START (3): one figure = two taus + Dr. Chen only
% plot(2*pi./omega, eta0, 'k', 'LineWidth', 1.5); %%% this is T which is T= 2*pi/omega plotting against T
plot(omega.^2 * RE / g, eta0, 'k', 'LineWidth', 1.5); %%% this is T which is T= 2*pi/omega plotting against T
hold on; % Add extra plots on same figure
% % % %
% % % % % Now plot the 3 CSV experimental points
% scatter(data_chen(:,1), data_chen(:,2), 30, 'r', 'o', 'filled');
% % scatter(data_liu(:,1), data_liu(:,2), 30, 'g', 's', 'filled');
% scatter(data_exp(:,1), data_exp(:,2), 30, 'b', '^', 'filled');
xlabel('$T$', 'Interpreter', 'latex');
ylabel('$|\eta / (iA)|$', 'Interpreter', 'latex');
title('Wave motion amplitude for $R_E = 138$', 'Interpreter', 'latex');
legend({'$\tau=0.2$','Model test'}, ...
'Interpreter','latex','Location','northwest');
grid on;
% xlim([5 35]); % Match the reference plot
% ylim([0 4.5]); % Optional: based on expected peak height
xlim([0 4]); % Match the reference plot
ylim([0 7]); % Optional: based on expected peak height
% === plotting section ===
% diary off
elapsedTime = toc;
disp(['Time consuming = ', num2str(elapsedTime), ' s']);
function out = dbesselk(l, z)
%DBESSELK Derivative of the modified Bessel function of the second kind
% out = dbesselk(l, z)
% Returns d/dz [K_l(z)] using the recurrence formula:
% K_l'(z) = -0.5 * (K_{l-1}(z) + K_{l+1}(z))
out = -0.5 * (besselk(l-1, z) + besselk(l+1, z));
end
function out = dbesselj(l, z)
%DBESSELJ Derivative of the Bessel function of the first kind
% out = dbesselj(l, z)
% Returns d/dz [J_l(z)] using the recurrence formula:
% J_l'(z) = 0.5 * (J_{l-1}(z) - J_{l+1}(z))
out = 0.5 * (besselj(l-1, z) - besselj(l+1, z));
end
function out = dbesseli(l, z)
%DBESSELI Derivative of the modified Bessel function of the first kind
% out = dbesseli(l, z)
% Returns d/dz [I_l(z)] using the recurrence formula:
% I_l'(z) = 0.5 * (I_{l-1}(z) + I_{l+1}(z))
out = 0.5 * (besseli(l-1, z) + besseli(l+1, z));
end
function out = dbesselh(l, z)
%DBESSELH Derivative of the Hankel function of the first kind
% out = dbesselh(l, z)
% Returns d/dz [H_l^{(1)}(z)] using the recurrence formula:
% H_l^{(1)'}(z) = 0.5 * (H_{l-1}^{(1)}(z) - H_{l+1}^{(1)}(z))
out = 0.5 * (besselh(l-1, 1, z) - besselh(l+1, 1, z));
end
function x = bessel0j(l,q,opt)
% a row vector of the first q roots of bessel function Jl(x), integer order.
% if opt = 'd', row vector of the first q roots of dJl(x)/dx, integer order.
% if opt is not provided, the default is zeros of Jl(x).
% all roots are positive, except when l=0,
% x=0 is included as a root of dJ0(x)/dx (standard convention).
%
% starting point for for zeros of Jl was borrowed from Cleve Moler,
% but the starting points for both Jl and Jl' can be found in
% Abramowitz and Stegun 9.5.12, 9.5.13.
%
% David Goodmanson
%
% x = bessel0j(l,q,opt)
k = 1:q;
if nargin==3 && opt=='d'
beta = (k + l/2 - 3/4)*pi;
mu = 4*l^2;
x = beta - (mu+3)./(8*beta) - 4*(7*mu^2+82*mu-9)./(3*(8*beta).^3);
for j=1:8
xnew = x - besseljd(l,x)./ ...
(besselj(l,x).*((l^2./x.^2)-1) -besseljd(l,x)./x);
x = xnew;
end
if l==0
x(1) = 0; % correct a small numerical difference from 0
end
else
beta = (k + l/2 - 1/4)*pi;
mu = 4*l^2;
x = beta - (mu-1)./(8*beta) - 4*(mu-1)*(7*mu-31)./(3*(8*beta).^3);
for j=1:8
xnew = x - besselj(l,x)./besseljd(l,x);
x = xnew;
end
end
end
% --- Local helper function for derivative of Bessel function ---
function dJ = besseljd(l, x)
dJ = 0.5 * (besselj(l - 1, x) - besselj(l + 1, x));
end
function v_D_val = v_D(N, Q, r, Z0d, wk, RI, l, d_vec, Znd, k, e_vec, chi)
% n = 0 mode
term1 = (d_vec(1) * Z0d * besselj(l, wk*r)) / (dbesselj(l, wk*RI));
% n >= 2 evanescent sum
sum1 = 0;
for nidx = 2:N
sum1 = sum1 + d_vec(nidx) * Znd(nidx) * besseli(l, k(nidx)*r) ...
/( dbesseli(l, k(nidx)*RI));
end
% q = 1..Q radial sum
sum2 = 0;
for qidx = 1:Q
sum2 = sum2 + e_vec(qidx) *chi(qidx)* besselj(l, chi(qidx)*r) / dbesselj(l, chi(qidx)*RI);
end
v_D_val = term1 + sum1 + sum2;
end
  10 commentaires
Javeria
Javeria le 8 Déc 2025
@Walter Roberson@Sam Chak what i have that we need to treat this non linear term iteratively to solve our system of equation;
The principle of iteration is the same. You may go in another way: 1) assuming your coefficient b1=0, you should have everything including v_D(r);
2) using a new a’( r)=a+b|v_D( r)| in place of a to solve all equations to get a new v_D( r); 3) to check the new v_D with previous one, if the difference is larger than your tolerance, continue the step 2).
David Goodmanson
David Goodmanson le 23 Déc 2025
Hi Sam, thanks for mentioning that I wrote a piece of the code, which I posted to Javeria on a previous thread. Javeria, if you see this, you can use the bessel0j code as you see fit but the nitpicky part of my nature wishes you had let the variabIe n (the bessel function order) alone, instead of changing it to l (small L). It's just a detail but in code generally, there are too many fonts where small L and capital i and the number 1 can be confused.

Connectez-vous pour commenter.

Réponses (4)

Torsten
Torsten le 9 Déc 2025
Modifié(e) : Torsten le 9 Déc 2025
Although it's slow, the full Newton method seems to work in this case (except for one value of X_kRE):
clear all;
close all;
clc;
tic;
% diary('iter_log.txt');
% diary on
% fprintf('Run started: %s\n', datestr(now));
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Parameters from your setup
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
old = false;
if old
N = 20; % Number of N modes
M = 20; % Number of M modes
Q = 20; % Number of Q modes
RI = 34.5; % Inner radius (m)
RE = 47.5; % Outer radius (m)
h = 200; % Water depth (m) % Interface depth (m)
d = 38; % Interface depth (m) % Draft (m)
tau = 0.1; % porosity ratio
X_kRE = linspace(0.0025, 0.16,100); %%%%%% this is now k_0
else
N = 10; % Number of N modes
M = 10; % Number of M modes
Q = 10; % Number of Q modes
RI = 0.4; % Inner radius (m)
RE = 0.8; % Outer radius (m)
h = 1.0; % Water depth (m) % Interface depth (m)
d = 0.3; % Interface depth (m) % Draft (m)
tau = 0.2; % porosity ratio
X_kRE = linspace(0.05, 4.5, 40); %%%%%% this is now k_0
end
b = h-d; % Distance from cylinder bottom to seabed (m)
g = 9.81; % Gravity (m/s^2)
gamma = 1.0; % discharge coefficient
b1 = (1 - tau) / (2 * gamma * tau^2); % nonlinear coeff
tol = 1e-4; % convergence tolerance
max_iter = 100; % max iterations
l = 0; % Azimuthal order
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% === Compute Bessel roots and scaled chi values ===
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
roots = bessel0j(l, Q); % Zeros of J_l(x)
chi = roots ./ RI; % chi_q^l = roots scaled by radius
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% ==============Initialize eta0 array before loop
%===============Initialize omega array before loop
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
eta0 = zeros(size(X_kRE)); % Preallocate
omega = zeros(size(X_kRE)); %%%% omega preallocate
% iters_used = zeros(length(X_kRE),1);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for idx = 1:length(X_kRE)
wk = X_kRE(idx); % dimensional wavenumber
omega(idx) = sqrt(g * wk * tanh(wk * h)); % dispersion relation
% fprintf('\n--- idx %3d (T=%6.3f s) ---\n', idx, 2*pi/omega(idx));
% drawnow;
%%%%%%%%%%%%% Compute group velocity %%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Cg(idx) = (g * tanh(wk * h) + g * wk * h * (sech(wk * h))^2) ...
* omega(idx) / (2 * g * wk * tanh(wk * h));
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% =======Compute a1 based on tau
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
a1(idx) = 0.93 * (1 - tau) * Cg(idx);
% dissip(idx) = a1(idx)/omega(idx);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%=============== Find derivative of Z0 at z = -d
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Z0d = (cosh(wk * h) * wk * sinh(wk * b)) / (2 * wk * h + sinh(2 * wk * h));
fun_alpha = @(x) omega(idx)^2 + g*x.*tan(x*h); %dispersion equation for evanescent modes
opts_lsq = optimset('Display','none'); % ← add this
%opts_lsq = [];
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,opts_lsq);
%%%%%%%%%%%%% Derivative of Z_n at z = -d %%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Znd(n) = -k(n) * (cos(k(n)*h) * sin(k(n)*b)) / (2*k(n)*h + sin(2*k(n)*h));
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
% % %%%%%%%%%% Write Matrix W %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
W = zeros(N,Q);
% %%W_{0q}
for q = 1:Q
r = exp(-2*chi(q)*b); % e^{-2*chi*b}
coth_chib = (1 + r) / (1 - r); % coth(chi*b) via decaying exp
W(1,q) = -(4*chi(q)/cosh(wk*h))*(( wk*sinh(wk*b)*coth_chib - chi(q)*cosh(wk*b) ) / ( wk^2 - chi(q)^2 ) ); %prof.chen
%
% % W(1, q) = - (4 * chi(q) / cosh(wk * h)) * (chi(q) * sinh(chi(q) * b) * cosh(wk * b) - wk * sinh(wk * b) * cosh(chi(q) * b)) ...
% % / ((chi(q)^2 - wk^2) * sinh(chi(q) * b));
end
% % %%% Write W_{nq}
for n = 2:N
for q = 1:Q
r = exp(-2*chi(q)*b); % e^{-2*chi*b}
coth_chib = (1 + r) / (1 - r); % coth(chi*b), stable
W(n,q) = -(4*chi(q)/cos(k(n)*h)) * ...
( ( k(n)*sin(k(n)*b)*coth_chib + chi(q)*cos(k(n)*b) ) ... %prof.chen
/ ( k(n)^2 + chi(q)^2 ) );
% % W(n, q) = - (4 * chi(q) / cos(k(n) * h)) * (chi(q) * sinh(chi(q) * b) * cos(k(n) * b) + k(n) * sin(k(n) * b) * cosh(chi(q) * b)) ...
% % / ((chi(q)^2 + k(n)^2) * sinh(chi(q) * b));
end
end
%%%%%%%%%%%%%%%%%%%%%%%Find E1_{q}%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
f2 = omega(idx)^2 / g;
E1 = zeros(Q,1);
S_q = zeros(Q,1);
for q = 1:Q
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% this expression is same as prof. Chen code
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
exp2chi_d = exp(-2 * chi(q) * d); % e^{-2*chi*d}
exp2chi_b = exp(-2 * chi(q) * b); % e^{-2*chi*(h-d)}
exp2chi_h = exp(-2 * chi(q) * h); % e^{-2*chi*h}
exp1chi_d = exp(-chi(q) * d); % e^{-chi*d}
% intermediate ratio (same order as Fortran)
BB = (f2 - chi(q)) * (exp2chi_b / exp2chi_d);
BB = BB + (f2 + chi(q)) * (1 + 2 * exp2chi_b);
BB = BB / (f2 - chi(q) + (f2 + chi(q)) * exp2chi_d);
BB = BB / (1 + exp2chi_b);
% final coefficient
E1(q) = 2 * (BB * exp2chi_d - 1);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% This expression is that one derived directly
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% % N0, D0, H0 exactly as in the math above
% N0 = (chi(q) - f2) + (chi(q) + f2) * exp2chi_h;
% D0 = (f2 - chi(q)) + (f2 + chi(q)) * exp2chi_d;
% H0 = 1 + exp2chi_b;
%
% % final coefficient (no add/subtract step)
% E1(q) = 2 * N0 / (D0 * H0);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% This one is the hyperbolic one form
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% num = chi(q)*cosh(chi(q)*h) - f2*sinh(chi(q)*h);
% den = (f2*cosh(chi(q)*d) - chi(q)*sinh(chi(q)*d)) * cosh(chi(q)*(h - d)); % or: * cosh(chi(q)*b)
% E1(q) = num / den;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%% This is the upper region velocity potential at z=0, write
%%%%%%%%%%%%%%%% from prof.chen code
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
CCmc = 1 + (exp2chi_b / exp2chi_d) + 2*exp2chi_b;
CCmc = CCmc / (1 + exp2chi_b);
CCmc = BB * (1 + exp2chi_d) - CCmc;
S_q(q) = CCmc * exp1chi_d; % this equals the summand value for this q
% --- direct hyperbolic (exact) ---
% % % C_hyp = ( chi(q)*cosh(chi(q)*h) - f2*sinh(chi(q)*h) ) ...
% % % / ( (f2*cosh(chi(q)*d) - chi(q)*sinh(chi(q)*d)) * cosh(chi(q)*b) );
% % %
% % % S_hyp(q) = C_hyp * cosh(chi(q)*d) + sinh(chi(q)*h) / cosh(chi(q)*b); % summand
end
% H(q): diagonal term (exponential form)
H1 = zeros(Q,1);
for q = 1:Q
% dissip =
r = exp(-2*chi(q)*b); % e^{-2*chi*b}
mid = 4*r/(1 - r^2) - E1(q); % 2/sinh(2*chi*b) - C_m
H1(q) = mid + 1i*a1(idx)*chi(q)/omega(idx); % mid + i*(a1/w)*chi
end
H = diag(H1);
%%%%%%%%%%matric G%%%%%%%%%%%%%%%%%%
G = zeros(Q, N); % Preallocate
for q = 1:Q
G(q,1) = 1i * 2*(a1(idx)/(omega(idx)*RI))* Z0d ...
* ( chi(q) / (chi(q)^2 - wk^2) )* ( besselj(l, wk*RI) / dbesselj(l, wk*RI) ); % this one is my calculated
end
% %G_qn
for q = 1:Q
for n = 2:N
G(q, n) = 1i * 2*(a1(idx)/(omega(idx)*RI))* Znd(n) ...
* ( chi(q) / (k(n)^2 + chi(q)^2) )* ( besseli(l, k(n)*RI) / dbesseli(l, k(n)*RI) ); % this one is my calculated
end
end
% % % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% %Define the right hand side vector
U = zeros(2*M + 2*N + Q, 1); % Full vector, M+N+M+N+Q elements
% % Block 1: U (size M = 4)
for i = 1:M
if i == 1
U(i, 1) = (sinh(wk*b))/((b*wk)*cosh(wk*h))*besselj(l, wk*RE) ; % Z_0^l
else
U(i, 1) = (2 *wk*(-1)^(i-1)*sinh(wk*b))/(b *(wk^2 +m(i)^2)*cosh(wk*h))*besselj(l, wk*RE); %Z_m^l
end
end
% % Block 2: Y (size N = 3)
for j = 1:N
if j == 1
U(j + M, 1) = -dbesselj(l, wk*RE) * (2*wk*h + sinh(2*wk*h)) /(cosh(wk*h)^2); % Y_0^l
else
U(j + M, 1) = 0; % Y_n^l
end
end
% % Block 3: X (size M = 4)
for i = 1:M
U(i + M + N, 1) = 0; % X_0^l, X_m^l
end
% Block 4: W (size N = 3)
for j = 1:N
U(j + M + N + M, 1) = 0; % W_0^l, W_n^l
end
for q = 1:Q
U(2*M + 2*N + q) = -(1i*b1/omega(idx)) * (2/(RI^2 * dbesselj(l, chi(q)*RI)));
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);
ZMQ = zeros(M,Q);
ZNQ = zeros(N,Q);
ZQM = zeros(Q,M);
ZQN = zeros(Q,N);
% Construct the full block matrix S
S = [ I_M, -A, ZMM, ZMN, ZMQ;
-B, I_N, -C, ZNN, ZNQ;
ZMM, ZMN, I_M, -D, ZMQ;
-E, ZNN, -F, I_N, -W;
ZQM, ZQN, ZQM, -G, H ];
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
solution_method = 0;
if solution_method == 0
if idx==1
psi0 = zeros(Q,1);
else
if exitflag == 1
psi0 = psi;
else
psi0 = zeros(Q,1);
end
end
f = @(psi)fun(psi,M,N,Q,S,U,Z0d,wk,RI,l,Znd,k,chi);
[psi,~,exitflag] = fsolve(f,psi0,optimset('Display','none'));
[~,d_vec,e_vec] = f(psi);
elseif solution_method == 1
converged = false;
psi = zeros(Q,1);
T_old = []; % <-- track the full solution from prev. iter
% T_old = zeros(2*M + 2*N + Q, 1); % previous unknowns (seed)
for iter = 1:max_iter
% (A) overwrite ONLY Block 5 of U from current psi
for q = 1:Q
U(2*M + 2*N + q) = -(1i*b1/omega(idx)) * (2/(RI^2 * dbesselj(l, chi(q)*RI))) * psi(q);
end
T = S \ U; % Solve the linear system
% T = pinv(S) * U; % Use pseudoinverse for stability
b_vec = T(1:M); % Coefficients b^l
a_vec = T(M+1 : M+N); % Coefficients a^l
c_vec = T(M+N+1 : 2*M+N); % Coefficients c^l
d_vec = T(2*M+N+1:2*M+2*N); % Coefficients d^l
e_vec = T(2*M+2*N+1:end); % (Q×1)
% (D) update psi for NEXT iteration from CURRENT coefficients
for q = 1:Q
integrand = @(r) abs(v_D(N,Q,r,Z0d,wk,RI,l,d_vec,Znd,k,e_vec,chi)) ...
.* v_D(N,Q,r,Z0d,wk,RI,l,d_vec,Znd,k,e_vec,chi) ...
.* besselj(l, chi(q)*r) .* r;
psi(q) = integral(integrand, 0, RI, 'AbsTol',1e-8, 'RelTol',1e-6);
end
% % % drawnow; % optional: flush output each iter
% === compact per-iteration print, only for T in [5,35] ===
% Tcur = 2*pi/omega(idx);
% if Tcur >= 5 && Tcur <= 35
% fprintf('idx %3d | iter %2d | ||psi|| = %.3e\n', idx, iter, norm(psi));
% drawnow;
% end
% % % % (E) convergence on ALL unknowns (T vs previous T)
% --- per-iteration diagnostics (optional) ---
if ~isempty(T_old)
diff_T = max(abs(T - T_old));
% Tcur = 2*pi/omega(idx);
% if Tcur >= 5 && Tcur <= 35
% fprintf('k0 idx %3d | iter %2d: max|ΔT| = %.3e, ||psi|| = %.3e\n', ...
% idx, iter, diff_T, norm(psi));
% drawnow; % ← add this to flush each iteration line too
% end
if diff_T < tol
converged = true;
break
end
end
T_old = T;
end % <<< end of iter loop
end
if solution_method==0
exitflag
elseif solution_method==1
converged
end
%T
% % ---- ONE summary line per frequency (place it here) ----
% iters_used(idx) = iter; % how many iterations this frequency needed
% fprintf('idx %3d (T=%6.3f s): iters=%2d\n', idx, 2*pi/omega(idx), iter);
% drawnow;
% -------------- end nonlinear iteration --------------
% % %%%%%%%%% 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
phiUell = 0;
for q =1:Q
phiUell = phiUell +e_vec(q)* S_q(q);
end
eta0(idx) = abs(term1+sum1+phiUell);
end
% % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% % % Plotting Data
% % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% %
% %
% figure(2); clf;
% >>> CHANGES START (3): one figure = two taus + Dr. Chen only
%if old
%plot(2*pi./omega, eta0, 'k', 'LineWidth', 1.5); %%% this is T which is T= 2*pi/omega plotting against T
%else
plot(omega.^2 * RE / g, eta0, 'k', 'LineWidth', 1.5); %%% this is T which is T= 2*pi/omega plotting against T
%end
hold on; % Add extra plots on same figure
% % % %
% % % % % Now plot the 3 CSV experimental points
% scatter(data_chen(:,1), data_chen(:,2), 30, 'r', 'o', 'filled');
% % scatter(data_liu(:,1), data_liu(:,2), 30, 'g', 's', 'filled');
% scatter(data_exp(:,1), data_exp(:,2), 30, 'b', '^', 'filled');
xlabel('$T$', 'Interpreter', 'latex');
ylabel('$|\eta / (iA)|$', 'Interpreter', 'latex');
title('Wave motion amplitude for $R_E = 138$', 'Interpreter', 'latex');
legend({'$\tau=0.2$','Model test'}, ...
'Interpreter','latex','Location','northwest');
grid on;
% xlim([5 35]); % Match the reference plot
% ylim([0 4.5]); % Optional: based on expected peak height
%xlim([0 4]); % Match the reference plot
% ylim([0 7]); % Optional: based on expected peak height
% === plotting section ===
% diary off
elapsedTime = toc;
disp(['Time consuming = ', num2str(elapsedTime), ' s']);
function [res,d_vec,e_vec] = fun(psi,M,N,Q,S,U,Z0d,wk,RI,l,Znd,k,chi)
for q = 1:Q
U(2*M + 2*N + q) = U(2*M + 2*N + q) * psi(q);
end
T = S\U;
d_vec = T(2*M+N+1:2*M+2*N); % Coefficients d^l
e_vec = T(2*M+2*N+1:end); % (Q×1)
for q = 1:Q
integrand = @(r) abs(v_D(N,Q,r,Z0d,wk,RI,l,d_vec,Znd,k,e_vec,chi)) ...
.* v_D(N,Q,r,Z0d,wk,RI,l,d_vec,Znd,k,e_vec,chi) ...
.* besselj(l, chi(q)*r) .* r;
res(q) = psi(q) - integral(integrand, 0, RI, 'AbsTol',1e-8, 'RelTol',1e-6);
end
end
function out = dbesselk(l, z)
%DBESSELK Derivative of the modified Bessel function of the second kind
% out = dbesselk(l, z)
% Returns d/dz [K_l(z)] using the recurrence formula:
% K_l'(z) = -0.5 * (K_{l-1}(z) + K_{l+1}(z))
out = -0.5 * (besselk(l-1, z) + besselk(l+1, z));
end
function out = dbesselj(l, z)
%DBESSELJ Derivative of the Bessel function of the first kind
% out = dbesselj(l, z)
% Returns d/dz [J_l(z)] using the recurrence formula:
% J_l'(z) = 0.5 * (J_{l-1}(z) - J_{l+1}(z))
out = 0.5 * (besselj(l-1, z) - besselj(l+1, z));
end
function out = dbesseli(l, z)
%DBESSELI Derivative of the modified Bessel function of the first kind
% out = dbesseli(l, z)
% Returns d/dz [I_l(z)] using the recurrence formula:
% I_l'(z) = 0.5 * (I_{l-1}(z) + I_{l+1}(z))
out = 0.5 * (besseli(l-1, z) + besseli(l+1, z));
end
function out = dbesselh(l, z)
%DBESSELH Derivative of the Hankel function of the first kind
% out = dbesselh(l, z)
% Returns d/dz [H_l^{(1)}(z)] using the recurrence formula:
% H_l^{(1)'}(z) = 0.5 * (H_{l-1}^{(1)}(z) - H_{l+1}^{(1)}(z))
out = 0.5 * (besselh(l-1, 1, z) - besselh(l+1, 1, z));
end
function x = bessel0j(l,q,opt)
% a row vector of the first q roots of bessel function Jl(x), integer order.
% if opt = 'd', row vector of the first q roots of dJl(x)/dx, integer order.
% if opt is not provided, the default is zeros of Jl(x).
% all roots are positive, except when l=0,
% x=0 is included as a root of dJ0(x)/dx (standard convention).
%
% starting point for for zeros of Jl was borrowed from Cleve Moler,
% but the starting points for both Jl and Jl' can be found in
% Abramowitz and Stegun 9.5.12, 9.5.13.
%
% David Goodmanson
%
% x = bessel0j(l,q,opt)
k = 1:q;
if nargin==3 && opt=='d'
beta = (k + l/2 - 3/4)*pi;
mu = 4*l^2;
x = beta - (mu+3)./(8*beta) - 4*(7*mu^2+82*mu-9)./(3*(8*beta).^3);
for j=1:8
xnew = x - besseljd(l,x)./ ...
(besselj(l,x).*((l^2./x.^2)-1) -besseljd(l,x)./x);
x = xnew;
end
if l==0
x(1) = 0; % correct a small numerical difference from 0
end
else
beta = (k + l/2 - 1/4)*pi;
mu = 4*l^2;
x = beta - (mu-1)./(8*beta) - 4*(mu-1)*(7*mu-31)./(3*(8*beta).^3);
for j=1:8
xnew = x - besselj(l,x)./besseljd(l,x);
x = xnew;
end
end
end
% --- Local helper function for derivative of Bessel function ---
function dJ = besseljd(l, x)
dJ = 0.5 * (besselj(l - 1, x) - besselj(l + 1, x));
end
function v_D_val = v_D(N, Q, r, Z0d, wk, RI, l, d_vec, Znd, k, e_vec, chi)
% n = 0 mode
term1 = (d_vec(1) * Z0d * besselj(l, wk*r)) / (dbesselj(l, wk*RI));
% n >= 2 evanescent sum
sum1 = 0;
for nidx = 2:N
sum1 = sum1 + d_vec(nidx) * Znd(nidx) * besseli(l, k(nidx)*r) ...
/( dbesseli(l, k(nidx)*RI));
end
% q = 1..Q radial sum
sum2 = 0;
for qidx = 1:Q
sum2 = sum2 + e_vec(qidx) *chi(qidx)* besselj(l, chi(qidx)*r) / dbesselj(l, chi(qidx)*RI);
end
v_D_val = term1 + sum1 + sum2;
end
  38 commentaires
Javeria
Javeria le 25 Déc 2025
@Torstenlicense('test','Distrib_Computing_Toolbox')
ans =
1
i tested my MATLAB
Torsten
Torsten le 25 Déc 2025
Modifié(e) : Torsten le 25 Déc 2025
First test whether the code works if you don't use the solution of the last step for a restart without parallel computing:
Replace
% Compute starting solution
x = X_kRE_min;
[converged,iter,resnorm,exitflag,sol,t,eta0,...
A_vec,B_vec,C_vec,D_vec,E_vec,...
AmA_vec,BmB_vec,FmC_vec,DmD_vec,WmE_vec] = ...
get_solution(N,M,Q,RI,RE,h,d,tau,x,g,gamma,l,...
solution_method,...
solver_method,...
TolX,TolF,itermax,...
TolX_nleq,TolFun_nleq,MaxIter_nleq,MaxFunEvals_nleq,...
AutoScaling_nleq,ComplexEqn_nleq,Updating_nleq,...
TypicalX_nleq,Jacobian_nleq,...
urelax,tol,max_iter,...
integral_method,RelTol_integral,AbsTol_integral,nr,...
solini);
if exitflag == 1 || converged == true
start = true;
else
start = false;
end
% If start in X_kRE_min is successful, continue solution for complete X_kRE vector
if start
disp("Starting computations ...");
disp(" ");
T = zeros(1,nXsteps+1);
Eta0 = zeros(1,nXsteps+1);
T(1) = t;
Eta0(1) = eta0;
Tfull(1) = t;
Eta0full(1)= eta0;
ns = 1;
%dx = stepX;
%Dx = stepX;
for i = 2:nXsteps+1
Xstart = X_kRE(i-1);
Xend = X_kRE(i);
%x = Xstart + Dx;
x = Xend;
dx = stepX;
sol0 = sol;
finished = false;
nsucc = 0;
while ~finished && dx >= stepX_min
[converged,iter,resnorm,exitflag,sol,t,eta0,...
A_vec,B_vec,C_vec,D_vec,E_vec,...
AmA_vec,BmB_vec,FmC_vec,DmD_vec,WmE_vec] = ... ...
get_solution(N,M,Q,RI,RE,h,d,tau,x,g,gamma,l,...
solution_method,...
solver_method,...
TolX,TolF,itermax,...
TolX_nleq,TolFun_nleq,MaxIter_nleq,MaxFunEvals_nleq,...
AutoScaling_nleq,ComplexEqn_nleq,Updating_nleq,...
TypicalX_nleq,Jacobian_nleq,...
urelax,tol,max_iter,...
integral_method,RelTol_integral,AbsTol_integral,nr,...
sol0);
if (solution_method <= 0 && ~converged) || ...
(solution_method > 0 && exitflag ~=1)
x
dx
eta0
if solution_method <= 0
converged
iter
else
exitflag
end
resnorm
disp(" ");
end
if exitflag == 1 || converged == true
ns = ns + 1;
Tfull(ns) = t;
Eta0full(ns) = eta0;
if abs(Xend - x) < 1e-8
finished = true;
break
end
%Dx = min(stepX,2*dx);
%Dx = dx;
nsucc = nsucc + 1;
if nsucc >= 3
dx = min(Xend-x,2*dx);
nsucc = 0;
else
dx = min(Xend-x,dx);
end
x = x + dx;
sol0 = sol;
else
nsucc = 0;
x = x - dx;
dx = dx/2;
%Dx = dx;
x = x + dx;
sol0 = sol0;
end
end
x
dx
eta0
disp(" ");
T(i) = t;
Eta0(i) = eta0;
end
% Write results to file ...
AA = [T.',Eta0.'];
% ... for prescribed grid points
fid = fopen('data_T', 'w+');
for i = 1:size(AA, 1)
fprintf(fid, '%f ', AA(i,:));
fprintf(fid, '\n');
end
fclose(fid);
AA = [Tfull.',Eta0full.'];
% ... for complete set of used grid points
fid = fopen('data_T_full', 'w+');
for i = 1:size(AA, 1)
fprintf(fid, '%f ', AA(i,:));
fprintf(fid, '\n');
end
fclose(fid);
hold on
plot(T,Eta0)
plot(Tfull,Eta0full)
hold off
grid on
size(Eta0)
size(Eta0full)
else
disp("Start was not successful.")
end
by
fid = fopen('data.txt', 'w+');
T = zeros(1,nXsteps+1);
Eta0 = zeros(1,nXsteps+1);
for i = 1:nXsteps+1
x = X_kRE(i);
[converged,iter,resnorm,exitflag,sol,t,eta0,...
A_vec,B_vec,C_vec,D_vec,E_vec,...
AmA_vec,BmB_vec,FmC_vec,DmD_vec,WmE_vec] = ...
get_solution(N,M,Q,RI,RE,h,d,tau,x,g,gamma,l,...
solution_method,...
solver_method,...
TolX,TolF,itermax,...
TolX_nleq,TolFun_nleq,MaxIter_nleq,MaxFunEvals_nleq,...
AutoScaling_nleq,ComplexEqn_nleq,Updating_nleq,...
TypicalX_nleq,Jacobian_nleq,...
urelax,tol,max_iter,...
integral_method,RelTol_integral,AbsTol_integral,nr,...
solini);
T(i) = t;
Eta0(i) = eta0;
fprintf(fid, '%f %f \n ',t,eta0);
end
fclose(fid);
plot(T,Eta0)
grid on
and see if it works for the case
N = 10; % Number of N modes
M = 10; % Number of M modes
Q = 10; % Number of Q modes
RI = 0.4; % Inner radius (m)
RE = 0.8; % Outer radius (m)
h = 1.0; % Water depth (m) % Interface depth (m)
d = 0.3; % Interface depth (m) % Draft (m)
tau = 0.2; % porosity ratio
g = 9.81; % Gravity (m/s^2)
gamma = 1.0; % discharge coefficient
l = 0; % Azimuthal order
X_kRE_min = 0.05
X_kRE_max = 4.5 %%%%%% this is now k_0 ;
nXsteps = 40
stepX = (X_kRE_max-X_kRE_min)/...
nXsteps;
X_kRE = linspace(X_kRE_min,X_kRE_max,nXsteps+1);
If yes, try parallelizing the loop with "parfor". I have no experience with parallel computing - so if you have problems, you will need to open a new question.

Connectez-vous pour commenter.


William Rose
William Rose le 31 Déc 2025
I am learning from your discussion. I also learned from the comments of @Sam Chak and @David Goodmanson. I like Sam's reminder to focus on the issues that most pertain to the original research question.
@Javeria, I strongly recommend that you accept @Torsten's answer.
  10 commentaires
Javeria
Javeria le 5 Jan 2026 à 14:10
@Torsten no i am not angry your suggestions are quite helpful for me and i am trying to get my required solution my capslock was on that time that why upper case letters was typed.
Torsten
Torsten le 7 Jan 2026 à 2:47
Modifié(e) : Torsten le 7 Jan 2026 à 14:50
I made some numerical experiments to find the reason why your fixed-point iteration works for one parameter set while it doesn't for the other.
Let's start simple to see the idea.
Let's compare the recursion x_(n+1) = 2*x_n with the recursion x_(n+1) = 0.5*x_n and an initial value x_0 ~=0. Then the first recursion diverges while the second converges to 0. Why is this so ?
A recursion is given by an iteration function f in the sense that x_(n+1) = f(x_n). For the above two cases, f(x) = 2*x for the first recursion and f(x) = 0.5*x for the second. A sufficient (and "almost" necessary) condition for convergence is |f'(x)| < 1 near the fixed point. This will ensure that if you start with x0 near the fixed point, the fixed point will "attract" the iterates.
If you generalize this fact to higher dimensions, the eigenvalues of the Jacobian matrix of f in the fixed point are the suitable indicator whether iteration will converge or not. So I computed the solution T* for the values of X_kRE with Newton iteration and deduced from it the Jacobian of the iteration function f at T*.
It turned out that in case of the parameter set where you got convergence, for all values of X_kRE, the Jacobian matrix of f in the T*'s only had eigenvalues of absolute values < 1 when you use the following setup:
M = 20; % Number of M modes
N = 20; % Number of N modes
Q = 20; % Number of Q modes
RI = 34.5; % Inner radius (m)
RE = 47.5; % Outer radius (m)
h = 200; % Water depth (m)
d = 38; % Interface depth (m) % Draft (m)
tau = 0.1; % porosity ratio
g = 9.81; % Gravity (m/s^2)
gamma = 1.0; % discharge coefficient
l = 0; % Azimuthal order
X_kRE_min = 0.0025;
X_kRE_max = 0.16; %%%%%% this is now k_0 ;
nXsteps = 100;
stepX = (X_kRE_max-X_kRE_min)/...
nXsteps;
X_kRE = linspace(X_kRE_min,X_kRE_max,nXsteps+1);
In case of the parameter set where you encounter divergence, the Jacobian at T* has eigenvalues in absolute value > 1 starting from X_kRE(15) when using the below setup (I didn't check when it ended) - and the fixed-point iteration started to diverge from X_kRE(17):
M = 20; % Number of M modes
N = 20; % Number of N modes
Q = 20; % Number of Q modes
RI = 0.345; % Inner radius (m)
RE = 0.475; % Outer radius (m)
h = 2.0; % Water depth (m) % Interface depth (m)
d = 0.38; % Interface depth (m) % Draft (m)
tau = 0.2; % porosity ratio
g = 9.81; % Gravity (m/s^2)
gamma = 1.0; % discharge coefficient
l = 0; % Azimuthal order
X_kRE_min = 0.1;%0.05;%2.0653;%0.05;
X_kRE_max = 4.2105;%20;%4.5;%2.5914;%4.5; %%%%%% this is now k_0 ;
nXsteps = 120;%40;%160;%40;
stepX = (X_kRE_max-X_kRE_min)/...
nXsteps;
X_kRE = linspace(X_kRE_min,X_kRE_max,nXsteps+1);
I hope I could convince you that using the fixed-point iteration approach independent of the parameter set doesn't make sense. It would be necessary to change the iteration function f in T_(n+1) = f(T_n) to have only eigenvalues of absolute value < 1 near the fixed points T*, but I have no idea how this could be achieved.

Connectez-vous pour commenter.


idris
idris le 9 Jan 2026 à 2:09
Modifié(e) : idris le 9 Jan 2026 à 2:16
Truly slow @Torsten, attached is the result
  6 commentaires
Javeria
Javeria le 12 Jan 2026 à 2:39
@idris no, i did not solve my problem
Javeria
Javeria le 12 Jan 2026 à 2:58
@Torsten see the below pdf and the equation number 24. They have the same non linear integral offcourse their expressions are changed but the non linear condition is same. They mentioned the method that how to solve it.

Connectez-vous pour commenter.


Javeria
Javeria le 12 Jan 2026 à 4:17
@Torsten i have the following code for the numerical integration but i am not sure whether it will work for my case and then how i can use it for my case.
clear
% Consider f(x) = cos(x) * exp(sin(x)), i.e. d/dx [exp(sin(x))]
f = @(x) cos(x) .* exp(sin(x));
a = 0.6; b = 2.2;
% The analytical evaluation is exp(sin(b)) - exp(sin(a))
int_analytical = exp(sin(b)) - exp(sin(a));
% Set the number of grid points
N = 101;
% Choose a method 'T', 'S', or 'G'. Need N odd for 'S', but N can be odd or
% even for 'T' and 'G'
method = 'T';
[x,w] = get_quadrature_weights(a, b, N, method);
int_numerical = sum(w.*f(x));
disp(['Analytical evaluation: ' num2str(int_analytical, 10)])
Analytical evaluation: 0.4857117348
disp(['Numerical evaluation: ' num2str(int_numerical, 10)])
Numerical evaluation: 0.4856852319
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [x,w] = get_quadrature_weights(a, b, N, method)
% get quadrature nodes (x) and weights (w) with N points over the
% interval a < x < b. The method can be 'T' (trapezium rule), 'S'
% (Simpson's rule) or 'G' (Gauss-Legendre quadrature).
% Trapezium rule: the spacing between points is equal
% Simpson's rule: N needs to be odd
if strcmp(method, 'T')
dx = (b-a)/(N-1);
x = linspace(a,b,N);
w = dx/2 * [1, 2*ones(1,N-2), 1];
elseif strcmp(method, 'S')
if mod(N,2) == 0
error('N must be odd to use Simpson''s rule')
return
end
dx = (b-a)/(N-1);
x = linspace(a,b,N);
n = (N-1)/2;
A = [4*ones(1,n-1); 2*ones(1,n-1)];
w = dx/3 * [1, A(:).', 4, 1];
elseif strcmp(method, 'G')
[x,w] = lgwt(N,a,b);
else
error('Method needs to be "T", "S" or "G"')
x = NaN; w = NaN;
end
end
% ------------------------------------------------------------------
function [x,w]=lgwt(N,a,b)
% lgwt.m
%
% This script is for computing definite integrals using Legendre-Gauss
% Quadrature. Computes the Legendre-Gauss nodes and weights on an interval
% [a,b] with truncation order N
%
% Suppose you have a continuous function f(x) which is defined on [a,b]
% which you can evaluate at any x in [a,b]. Simply evaluate it at all of
% the values contained in the x vector to obtain a vector f. Then compute
% the definite integral using sum(f.*w);
%
% Written by Greg von Winckel - 02/25/2004
N=N-1;
N1=N+1; N2=N+2;
xu=linspace(-1,1,N1)';
% Initial guess
y=cos((2*(0:N)'+1)*pi/(2*N+2))+(0.27/N1)*sin(pi*xu*N/N2);
% Legendre-Gauss Vandermonde Matrix
L=zeros(N1,N2);
% Derivative of LGVM
Lp=zeros(N1,N2);
% Compute the zeros of the N+1 Legendre Polynomial
% using the recursion relation and the Newton-Raphson method
y0=2;
% Iterate until new points are uniformly within epsilon of old points
while max(abs(y-y0))>eps
L(:,1)=1;
Lp(:,1)=0;
L(:,2)=y;
Lp(:,2)=1;
for k=2:N1
L(:,k+1)=( (2*k-1)*y.*L(:,k)-(k-1)*L(:,k-1) )/k;
end
Lp=(N2)*( L(:,N1)-y.*L(:,N2) )./(1-y.^2);
y0=y;
y=y0-L(:,N2)./Lp;
end
% Linear map from[-1,1] to [a,b]
x=(a*(1-y)+b*(1+y))/2;
% Compute the weights
w=(b-a)./((1-y.^2).*Lp.^2)*(N2/N1)^2;
end
  32 commentaires
Torsten
Torsten le 28 Jan 2026 à 13:36
Modifié(e) : Torsten le 28 Jan 2026 à 19:26
Even in the code I provided, I can't get convergence for all values of T immediately. And your expectations are too high if you think there could be a trick that could make this possible. But until now, I succeed to get a converged solution in all requested points by inserting one or multiple T values in between T(i-1) (where I got a converged solution) and T(i) (where I didn't get a converged solution) always starting with the last converged solution for the unknowns as initial guess. Moreover, I applied the real Newton's method for real(f) and imag(f) instead of the complex Newton's method to reduce the number of necessary iterations and I solved for T(2*M+2*N+1:2*M+2*N+Q) instead of psi(1:Q). But as you can see from my code: these modifications need some programming effort - it's not a one-liner that could be added in your code from above.
Enclosed are the results for your last computational case with H = 10 and H = 20 using MATLAB's "integral" as integration method for comparison using the code I provided.
M = 20; % Number of M modes
N = 20; % Number of N modes
Q = 20; % Number of Q modes
RI = 34.5; % Inner radius (m)
ratio = 47.5/34.5; % Ratio R_E / R_I
RE = ratio*RI; % outer radius (m)
h = 200/34.5 * RI; % Water depth (m)
d = 38/34.5 * RI; % Interface depth (m)
tau = 0.2; % porosity ratio
g = 9.81; % Gravity (m/s^2)
gamma = 1.0; % discharge coefficient
l = 0; % Azimuthal order
nXsteps = 199;
T1 = linspace(10, 19, nXsteps+1); % Define wave period range [s]
% T1 = linspace(0.5, 3.5, 50); % Scaled period range [s] [5/10 to 35/10]
omega = 2 * pi ./ T1; % Angular frequency [rad/s]
X_kRE = zeros(nXsteps+1,1); % Preallocate k_0 array
for i = 1:nXsteps+1
fun = @(x) omega(i)^2 - g * x * tanh(x * h);
X_kRE(i) = fsolve(fun,1,optimset('Display','none'));
end
A10 = dlmread('data_H=10_integral_focussed.txt');
figure(1)
plot(A10(:,1),A10(:,2))
xlabel('T')
ylabel('eta0')
title('H = 10')
grid on
A20 = dlmread('data_H=20_integral_focussed.txt');
figure(2)
plot(A20(:,1),A20(:,2))
xlabel('T')
ylabel('eta0')
title('H = 20')
grid on
Torsten
Torsten le 29 Jan 2026 à 13:50
That's why in the code I provided, I use the psi I get in iteration (i-1) as start point in iteration i.
As far as I can see, you don't do that. You start with psi=zeros(Q,1) for each value of T:
% --- Newton parameters
psi = zeros(Q,1); % initial guess (linear solution)
tol_NR = 1e-6;
max_NR = 100;
nr = 160;
R = linspace(0,RI,nr);
I meant you should try
% --- Newton parameters
if idx == 1 || it == max_NR
psi = zeros(Q,1); % initial guess (linear solution)
end
tol_NR = 1e-6;
max_NR = 100;
nr = 160;
R = linspace(0,RI,nr);
to have a better initial guess from the psi-vector of the last solution for T(idx-1) - thus reducing the number of iterations needed for T(idx).

Connectez-vous pour commenter.

Catégories

En savoir plus sur Linear Algebra 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!

Translated by