numerical Instabilities for bessel functions
Ältere Kommentare anzeigen
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 Kommentare
Javeria
am 8 Dez. 2025
Sam Chak
am 8 Dez. 2025
Hi @Javeria
Your code is quite lengthy. If it takes more than a minute to run on this platform, you are likely to receive the "unauthorized to perform this operation" message. Have you identified the specific portion of the Bessel function in your code that incurs numerical instabilities? I also notice that part of your code incorporates contributions from @David Goodmanson.
Javeria
am 8 Dez. 2025
Javeria
am 8 Dez. 2025
Walter Roberson
am 8 Dez. 2025
I am not convinced that the problem is due to instability in bessel functions.
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;
after awhile, d_vec and e_vec both grow to roughly +/- 1e172 +/- 1e172i . As a result, v_D returns results in roughly the same range. Multiply two such results together and you get overflow to +/- infinity. The besselj values remain modest.
Javeria
am 8 Dez. 2025
Sam Chak
am 8 Dez. 2025
Hi @Javeria
@Walter Roberson noted that both d_vec and e_vec grow to extremely large values and that they are related to the linear system S*T = U. Is the system mathematically guaranteed to have solutions? If so, the very large values in the system are probably causing numerical overflow.
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)
Javeria
am 8 Dez. 2025
David Goodmanson
am 23 Dez. 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.
Antworten (4)
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 Kommentare
Javeria
am 9 Dez. 2025
Javeria
am 9 Dez. 2025
Javeria
am 9 Dez. 2025
Javeria
am 9 Dez. 2025
Your iteration is some kind of fixed-point iteration, Newton's method uses derivative information to get the roots of the nonlinear system. So, Newton's method uses more information about the functions in its iteration process and thus usually has a larger region of convergence than the fixed-point iteration. Further, "fsolve" (the nonlinear solver I used) uses a damped Newton method that is even more sophisticated than the classical Newton method. I think it will lead you nowhere to speculate about a particular reason for divergence in your special case.
And you shouldn't interprete too much into the "second peak" in the result's graph. If you look at the exitflag from the nonlinear solver, it changes here from 1 to 3 (I think) in the point before the second peak (and has a lot of trouble converging in the points nearby as you can see from the long iteration times for these points). So I guess the result before the second peak (the T ~ 0.8 result) is not trustworthy, and at the moment I don't know how to make the solver converge for this particular value of X_kRE.
Because I start the iteration in point i from a trustworthy result in step (i-1), I tried with X_kRE = linspace(0.05, 4.5, 80); instead of X_kRE = linspace(0.05, 4.5, 40); to be nearer to the new solution when restarting, but this even gives exitflag = -3 for two points in the critical region. For the meaning of the exitflags, see
Javeria
am 9 Dez. 2025
You get a converged solution for both cases - with the exception of one input point - if you use a sophisticated root finder as in the code above.
I tried to explain that fixed-point methods usually have a smaller radius of convergence around the solution than Newton-like methods.
Since both methods (your method and Newton's method) start from the zero solution and your fixed-point method fails while Newton's method succeeds, this is the most probable reason for the different behaviour.
The physics for the intermediate water depth case seems to be numerically more challenging than the deep water case and needs a more sophisticated solver for the nonlinear contribution. Maybe you find a reason in the physics of the two problems why this is the case. If you look at the problem only from the numerical standpoint, it's too complex to identify the reason for failure with your solution method (at least for me).
Numerics is trial and error: if you know that your equations are correct and your solution method fails, you usually try another solution method. If this works, you don't ask why the first failed because it usually cannot be answered easily.
By the way: Newton's method also works for the deep water conditions (although it's again slower than fixed-point iteration which also succeeds in this case).
Javeria
am 10 Dez. 2025
Sam Chak
am 10 Dez. 2025
Hi @Javeria
I couldn't help because as I haven't followed your previous discussions. However, if your studies focus on wave motion, your aim should be to select the most effective and efficient solver. If your research involves successfully applying a novel numerical method to solve wave motion, you need not abandon the proposed method simply because it is challenging to solve the system at intermediate water depths.
In fact, Newtonian mechanics, which assumes that an object's mass remains constant regardless of its velocity (𝑣), also fails at speeds approaching the speed of light (𝑐). In reality, mass increases as an object approaches the speed of light, approaching infinity as 𝑣 nears 𝑐.
Moreover, your fixed-point iteration method outperforms Newton's method in the deep water case. Therefore, you may consider a hybrid solver approach that mimics Continuously Variable Transmission (CVT) smoothness in automatically selecting Newton's method for intermediate water depths and the fixed-point iteration method for the deep water case.
Javeria
am 10 Dez. 2025
I don't think that changing the code to compute the integral
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);
can help to solve the present problem.
What you have to do is to solve a system of equations
S*T = U(T)
for T.
What you do in your code is to perform a fixed-point iteration
T_(n+1) = S\U(T_n)
and hope that T_n will converge to some vector T.
The equivalent formulation would be to use Newton's method to solve
S*T - U(T) = 0
for T (I do something slightly different since I don't solve for T, but for psi).
Javeria
am 11 Dez. 2025
Torsten
am 11 Dez. 2025
Out of interest: how did you change your code to get the iteration converge ?
Javeria
am 11 Dez. 2025
David Goodmanson
am 11 Dez. 2025
Hi Sam,
The concept of equivalent relativistic mass is ok for colloquial arguments, but physics has kind of moved away from that concept. I think that's mostly because of the inclusion of general relativity into the mainstream, where people think about spacetime and other funcamental quantities like (E,p) that transform relativistically. If m is the invariant rest mass, then
E = mc^2/sqrt(1-v^2/c^2)
p = mv/sqrt(1-v^2/c^2)
E^2 = p^2c^2 + m^2c^4
p = 0 ==> E = mc^2 % object at rest
It's algebraically correct but a bit dubious to conclude from the first and last lines that
E/c^2 = m/sqrt(1-v^2/c^2) = m_equiv.
m_equiv is a quantity that kind of sits on its own. Without it, the speed-up argument is the same but different: To take a massive (rest mass) object up to v=c in a given inertial frame you have to pump in an unllimited amount of energy.
Thsnks but at the end we have the converged solution for each X_kRE(this is actually k0 in my code).
Ok, I misunderstood this comment from your side. I interpreted it that you already solved your problem with the help of the attached article.
Which part of your code could be affected by the article ?
Javeria
am 12 Dez. 2025
I used the modified code to compute eta0 as a function of X_kRE with a finer resolution around the peak. I got the following result:
A = readmatrix('data.txt');
plot(A(:,1),A(:,2))
xlabel('$T$', 'Interpreter', 'latex');
ylabel('$|\eta / (iA)|$', 'Interpreter', 'latex');
title('Wave motion amplitude for $R_E = 138$', 'Interpreter', 'latex');
legend({'$\tau=0.2$'}, ...
'Interpreter','latex','Location','northwest');
grid on
I left your computational part unchanged, but had to add some code in order to make it work according to my ideas.
In the additionally attached "program1.m", you can also use your fixed-point iteration if you set "solution_method = 0". This option will not work in "program.m". So use "program1.m" instead of "program.m" since it is an extension of "program.m".
If you have questions, you can of course come back and ask.
Javeria
am 18 Dez. 2025
Javeria
am 22 Dez. 2025
Javeria
am 22 Dez. 2025
If you change "ode=true" to "old=false" and leave "solution_method = 0", you try to run your fixed-point iteration method together with these parameters:
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
This will not work and that's why you came here to ask for a solution.
Set "solution_method" to 1 or 2 to make it work.
Concerning the integration method, I use "quadgk", a Gauss-Kronos quadrature rule, instead of "integral". In some cases, this worked better than "integral". If you want to return to "integral", comment out the lines where "quadgk" is called and replace them by the "integral" calls (that are still present at these code positions).
I additionally added a simple trapezoidal rule to evaluate the integral which can be chosen by using "method_integral = 1" and an appropriate number of discretization points "nr".
If you have further questions, please include the complete parameter setting part:
solution_method = 0; % 0: Fixed-point iteration
% 1: Newton method (psi-based)
% 2: Newton method (minimal-T-based)
% 3: Newton method (medium-T-based) (too slow, not recommended)
% 4: Newton method (full-T-based) (too slow, not recommended)
method_integral = 0; % 0: Use integral
% 1: Use trapezoidal rule
nr = 300; % if method_integral = 1: number of discretization points in r-direction
% Used if method_integral = 0
RelTol_integral = 1e-6; % Relative error tolerance in integral solving (default value: 1e-6)
AbsTol_integral = 1e-10; % Absolute error tolerance in integral solving (default value: 1e-10)
% Used if solution_method = 0
urelax = 0.01; % Underrelaxation factor
tol = 1e-4; % convergence tolerance
max_iter = 500; % max iterations
% Used if solution_method ~= 0
TolX_nleq = 1e-6; % Relative error tolerance in integral solving (default value: 1e-6)
TolFun_nleq = 1e-6; % Absolute error tolerance in integral solving (default value: 1e-6)
AutoScaling_nleq = 'off'; % default value: 'off'
ComplexEqn_nleq = 'off'; % default value: 'off'
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)
d = 38; % Interface depth (m) % Draft (m)
tau = 0.1; % porosity ratio
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);
stepX_min = stepX/2^6;
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_min = 2.0653;%0.05;%2.0653;
X_kRE_max = 2.5914;%4.5;%2.5914; %%%%%% this is now k_0 ;
nXsteps = 160;%40;%160;
stepX = (X_kRE_max-X_kRE_min)/...
nXsteps;
X_kRE = linspace(X_kRE_min,X_kRE_max,nXsteps+1);
stepX_min = stepX/2^15;
end
Usually, it should not be necessary to change parts of the code below these lines and so it will help to identify the problem.
Javeria
am 23 Dez. 2025
I developed the code under Octave, and the exitflags of "fsolve" under octave are different. So please only use the below code with solver_method = 11 or solver_method = 12 as these options invoked code that is "fsolve" - independent. But I strongly recommend to make the "fsolve" option work under MATLAB, too (the code below works without problems under Octave with solver_method = 21 and solver_method = 22):
Under MATLAB answers, the code below needed 32 seconds for 6 values of the X_kRE array. I requested 161 - thus the approximate runtime should be 32*161/6 seconds (maybe some more, because more iterations are needed for bigger values of X_kRE). On my Laptop, it took 3016 sec under octave. But you can see the progress on the terminal.
Please run the below code and report if it works.
After this, remove the lines
nXsteps = 5;
X_kRE=X_kRE(1:6);
to run the full second case with 160 instead of 6 output points in the critical region. How long does it take on your computer ?
main()
function main
tic
%warning('off')
solution_method = 2; % -4: Fixed-point iteration (full-T-based) (unstable, not recommended)
% -3: Fixed-point iteration (medium-T-based) (unstable, not recommended)
% -2: Fixed-point iteration (minimal-T-based) (unstable, not recommended)
% -1: Fixed-point iteration (psi-based) (unstable, not recommended)
% 0: Fixed-point iteration (Ali) (unstable, not recommended)
% 1: Newton method (psi-based)
% 2: Newton method (minimal-T-based)
% 3: Newton method (medium-T-based) (too slow, not recommended)
% 4: Newton method (full-T-based) (too slow, not recommended)
solver_method = 12; % if solution_method > 0:
% 11: Use self-programmed Newton method to solve the nonlinear
% complex-valued system of equations directly
% 12: Use self-programmed Newton method to solve the nonlinear
% complex-valued system of equations by splitting in real and imaginary part
% 21: Use "fsolve" to solve the nonlinear complex-valued system
% of equations directly
% 22: Use "fsolve" to solve the nonlinear complex-valued system
% of equations by splitting in real and imaginary part
% if solution_method < 0:
% 11 or 21 : Iterate the complex-valued solution vector
% 12 or 22 : Iterate the real and complex parts of the
% solution vector separately
% Used if solution_method > 0 & solver_method = 11 or 12:
TolX = 1e-6; % if solver_method = 2: Error tolerance in iterate
TolF = 1e-6; % if solver_method = 2: Error tolerance in function value
itermax = 50; % if solver_method = 2: Maximum number of Newton iterations
% Used if solution_method > 0 & solver_method = 21 or 22:
TolX_nleq = 1e-6; % Error tolerance in iterate (default value: 1e-6)
TolFun_nleq = 1e-6; % Error tolerance in function value (default value: 1e-6)
MaxIter_nleq = 400; % Maximum number of Newton iterations (default value: 400)
MaxFunEvals_nleq = []; % Maximum number of function evaluations (default value: [])
AutoScaling_nleq = 'off'; % default value: 'off'
ComplexEqn_nleq = 'off'; % default value: 'off'
Updating_nleq = 'off'; % Use Broyden updates to approximate the Jacobian (default value: 'off')
TypicalX_nleq = []; % Typical order of the solution (default value : [])
Jacobian_nleq = 'off'; % User-supplied Jacobian (default value: 'off')
% Used if solution_method = 0
urelax = 1; % Underrelaxation factor
tol = 1e-6; % convergence tolerance
max_iter = 50; % max iterations
% Set integration method
integral_method = 0; % 0: Use integral, 1: Use trapezoidal rule
RelTol_integral = 1e-6; % if integral_method = 0: Relative error tolerance in integral solving (default value: 1e-6)
AbsTol_integral = 1e-10; % if integral_method = 0: Absolute error tolerance in integral solving (default value: 1e-10)
nr = 300; % if integral_method = 1: number of discretization points in r-direction
% Set computational case
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)
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);
stepX_min = stepX/2^6;
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
g = 9.81; % Gravity (m/s^2)
gamma = 1.0; % discharge coefficient
l = 0; % Azimuthal order
X_kRE_min = 2.0653;%0.05
X_kRE_max = 2.5914;%4.5 %%%%%% this is now k_0 ;
nXsteps = 160;%40
stepX = (X_kRE_max-X_kRE_min)/...
nXsteps;
X_kRE = linspace(X_kRE_min,X_kRE_max,nXsteps+1);
nXsteps = 5;
X_kRE=X_kRE(1:6);
stepX_min = stepX/2^15;
end
% Initialize solution vectors
if solver_method == 11 || solver_method == 21
if abs(solution_method) <= 2
solini = zeros(Q,1);
elseif abs(solution_method) == 3
solini = zeros(N+Q,1);
elseif abs(solution_method) == 4
solini = zeros(2*N+2*M+Q,1);
end
elseif solver_method == 12 || solver_method == 22
if abs(solution_method) <= 2
solini = zeros(2*Q,1);
elseif abs(solution_method) == 3
solini = zeros(2*(N+Q),1);
elseif abs(solution_method) == 4
solini = zeros(2*(2*N+2*M+Q),1);
end
end
% 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
toc
end
function [converged,iter,resnorm,exitflag,sol,xeta0,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_kRE,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);
converged = false;
iter = 0;
resnorm = 0;
exitflag = 0;
%clear all;
%close all;
%clc;
%tic;
% diary('iter_log.txt');
% diary on
% fprintf('Run started: %s\n', datestr(now));
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Parameters from your setup
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%if solution_method == 1 || solution_method == 2 || ...
% solution_method == 3 || solution_method == 4
% first = true;
%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
%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);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
A_vec = zeros(N,length(X_kRE));
B_vec = zeros(M,length(X_kRE));
C_vec = zeros(M,length(X_kRE));
D_vec = zeros(N,length(X_kRE));
E_vec = zeros(Q,length(X_kRE));
AmA_vec = zeros(M,length(X_kRE));
BmB_vec = zeros(N,length(X_kRE));
FmC_vec = zeros(N,length(X_kRE));
DmD_vec = zeros(M,length(X_kRE));
WmE_vec = zeros(N,length(X_kRE));
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
% % 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 ];
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if solution_method < 0 % Fixed-point iterations
Y = sol0;
f = @(Y) wrapper(Y,S,U,M,N,Q,Z0d,wk,RI,l,Znd,k,chi,b1,omega,idx,...
solution_method,solver_method,...
integral_method,RelTol_integral,AbsTol_integral,nr);
for iter = 1:max_iter
[Y_new,a_vec,b_vec,c_vec,d_vec,e_vec] = f(Y);
diff_Y = max(abs(Y_new - Y));
if diff_Y < tol
converged = true;
break
end
Y = Y_new;
end
sol = Y_new;
%
%Jac = JacMatrix(f,sol);
%rho = norm(Jac)
resnorm = diff_Y;
elseif solution_method == 0 % Ali
%converged = false;
%psi = zeros(Q,1);
psi = sol0;
%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;
if integral_method == 0
psi(q) = quadgk(integrand, 0, RI, ...
'RelTol',RelTol_integral,'AbsTol',AbsTol_integral);
% psi(q) = integral(integrand, 0, RI, ...
% 'RelTol',RelTol_integral,'AbsTol',AbsTol_integral);
elseif integral_method == 1
R = linspace(0,RI,nr);
psi(q) = trapz(R,integrand(R));
end
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
T_old = urelax*T + (1-urelax)*T_old;
end
end % <<< end of iter loop
sol = psi;
resnorm = diff_T;
elseif solution_method > 0 % Newton-based formulations
if solver_method == 11 || solver_method == 12
f = @(Y) wrapper(Y,S,U,M,N,Q,Z0d,wk,RI,l,Znd,k,chi,b1,omega,idx,...
solution_method,solver_method,...
integral_method,RelTol_integral,AbsTol_integral,nr);
df = @(Y) JacMatrix(@(Z)f(Z),Y);
Yold = sol0;
iter = 0;
while iter < itermax
iter = iter + 1;
[dfY,fY,a_vec,b_vec,c_vec,d_vec,e_vec] = df(Yold);
dY = -dfY\fY;
Ynew = Yold + dY;
ErrorX = max(abs(dY));
ErrorF = max(abs(fY));
if ErrorX < TolX && ErrorF < TolF
exitflag = 1;
break
end
Yold = Ynew;
end
iter
sol = Ynew;
resnorm = ErrorF;
elseif solver_method == 21 || solver_method == 22
f = @(Y) wrapper(Y,S,U,M,N,Q,Z0d,wk,RI,l,Znd,k,chi,b1,omega,idx,...
solution_method,solver_method,...
integral_method,RelTol_integral,AbsTol_integral,nr);
df = [];
if strcmp(Jacobian_nleq,'on')
df = @(Y)JacMatrix(@(Z)f(Z),Y);
% df = @(Y)jacobianest(@(Z)f(Z),Y);
% df = @(Y)numgradient('wrapper',{Y,S,U,M,N,Q,Z0d,wk,RI,l,Znd,k,chi,b1,omega,idx,...
% solution_method,...
% integral_method,RelTol_integral,AbsTol_integral,nr});
end
options = optimset('Display','none','TolX',TolX_nleq,'TolFun',TolFun_nleq,...
'MaxIter',MaxIter_nleq,'MaxFunEvals',MaxFunEvals_nleq,...
'AutoScaling',AutoScaling_nleq,'ComplexEqn',ComplexEqn_nleq,...
'Updating',Updating_nleq,'TypicalX',TypicalX_nleq,...
'Jacobian',Jacobian_nleq);
[sol,resnorm,exitflag,~,Jac_fsolve] = fsolve({f,df},sol0,options);
[~,a_vec,b_vec,c_vec,d_vec,e_vec] = f(sol);
resnorm = norm(resnorm);
%Jac_fsolve
%Jac = JacMatrix(@(Z)f(Z),sol);
%Jac = jacobianest(@(Z)f(Z),sol)
%Jac = numgradient('wrapper',{sol,S,U,M,N,Q,Z0d,wk,RI,l,Znd,k,chi,b1,omega,idx,...
% solution_method,solver_method,...
% integral_method,RelTol_integral,AbsTol_integral,nr})
end
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);
xeta0(idx) = omega(idx)^2*RE/g;
A_vec(:,idx) = a_vec;
B_vec(:,idx) = b_vec;
C_vec(:,idx) = c_vec;
D_vec(:,idx) = d_vec;
E_vec(:,idx) = e_vec;
AmA_vec(:,idx) = A*a_vec;
BmB_vec(:,idx) = B*b_vec;
FmC_vec(:,idx) = F*c_vec;
DmD_vec(:,idx) = D*d_vec;
WmE_vec(:,idx) = W*e_vec;
%if solution_method == 0
% idx
% converged
% iter
% if ~converged
% eta0(idx) = NaN;
% endif
%elseif solution_method == 1 || solution_method == 2 || ...
% solution_method == 3 || solution_method == 4
% idx
% if exitflag == 1
% exitflag,exitflag_check
% norm(resnorm),norm(resnorm_check)
% else
% exitflag
% norm(resnorm)
% endif
% if exitflag ~= 1
% eta0(idx) = NaN;
% endif
end
%eta0(idx)
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']);
%end
function [fY,a_vec,b_vec,c_vec,d_vec,e_vec] = ...
wrapper(Y,S,U,M,N,Q,Z0d,wk,RI,l,Znd,k,chi,b1,omega,idx,...
solution_method,solver_method,...
integral_method,RelTol_integral,AbsTol_integral,nr)
if solver_method == 11 || solver_method == 21
[fY,a_vec,b_vec,c_vec,d_vec,e_vec] = ...
compute_fY(Y,S,U,M,N,Q,Z0d,wk,RI,l,Znd,k,chi,b1,omega,idx,...
solution_method,...
integral_method,RelTol_integral,AbsTol_integral,nr);
elseif solver_method == 12 || solver_method == 22
Yc = Y(1:end/2) + 1i*Y(end/2+1:end);
[fYc,a_vec,b_vec,c_vec,d_vec,e_vec] = ...
compute_fY(Yc,S,U,M,N,Q,Z0d,wk,RI,l,Znd,k,chi,b1,omega,idx,...
solution_method,...
integral_method,RelTol_integral,AbsTol_integral,nr);
fY = [real(fYc);imag(fYc)];
end
end
function [fY,a_vec,b_vec,c_vec,d_vec,e_vec] = ...
compute_fY(Y,S,U,M,N,Q,Z0d,wk,RI,l,Znd,k,chi,b1,omega,idx,...
solution_method,...
integral_method,RelTol_integral,AbsTol_integral,nr)
if solution_method == -1 || solution_method == 1
psi = Y;
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;
a_vec = T(M+1:M+N);
b_vec = T(1:M);
c_vec = T(M+N+1:2*M+N);
d_vec = T(2*M+N+1:2*M+2*N);
e_vec = T(2*M+2*N+1:2*M+2*N+Q);
fpsi = zeros(size(psi));
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;
if integral_method == 0
fpsi(q) = quadgk(integrand, 0, RI, ...
'RelTol',RelTol_integral,'AbsTol',AbsTol_integral);
% fpsi(q) = integral(integrand, 0, RI, ...
% 'RelTol',RelTol_integral,'AbsTol',AbsTol_integral);
elseif integral_method == 1
R = linspace(0,RI,nr);
fpsi(q) = trapz(R,integrand(R));
end
end
if solution_method == -1
fpsi = fpsi;
elseif solution_method == 1
fpsi = psi - fpsi;
end
fY = fpsi;
elseif solution_method == -2 || solution_method == 2
T = Y;
Ur = U(1:2*M+2*N) - S(1:2*M+2*N,2*M+2*N+1:2*M+2*N+Q)*T;
Tr = S(1:2*M+2*N,1:2*M+2*N)\Ur;
a_vec = Tr(M+1:M+N);
b_vec = Tr(1:M);
c_vec = Tr(M+N+1:2*M+N);
d_vec = Tr(2*M+N+1:2*M+2*N);
e_vec = T;
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;
if integral_method == 0
psi = quadgk(integrand, 0, RI, ...
'RelTol',RelTol_integral,'AbsTol',AbsTol_integral);
% psi = integral(integrand, 0, RI, ...
% 'RelTol',RelTol_integral,'AbsTol',AbsTol_integral);
elseif integral_method == 1
R = linspace(0,RI,nr);
psi = trapz(R,integrand(R));
end
U(2*M+2*N+q) = -(1i*b1/omega(idx)) * (2/(RI^2 * dbesselj(l, chi(q)*RI))) * psi;
end
if solution_method == -2
Unr = U(2*M+2*N+1:2*M+2*N+Q) - S(2*M+2*N+1:2*M+2*N+Q,1:2*M+2*N)*Tr;
Tnr = S(2*M+2*N+1:2*M+2*N+Q,2*M+2*N+1:2*M+2*N+Q)\Unr;
fT = Tnr;
elseif solution_method == 2
fT = S(2*M+2*N+1:2*M+2*N+Q,1:2*M+2*N)*Tr(1:2*M+2*N) + ...
S(2*M+2*N+1:2*M+2*N+Q,2*M+2*N+1:2*M+2*N+Q)*T - ...
U(2*M+2*N+1:2*M+2*N+Q);
end
fY = fT;
elseif solution_method == -3 || solution_method == 3
T = Y;
Ur = U(1:2*M+N) - S(1:2*M+N,2*M+N+1:2*M+2*N+Q)*T;
Tr = S(1:2*M+N,1:2*M+N)\Ur;
a_vec = Tr(M+1:M+N);
b_vec = Tr(1:M);
c_vec = Tr(M+N+1:2*M+N);
d_vec = T(1:N);
e_vec = T(N+1:N+Q);
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;
if integral_method == 0
psi = quadgk(integrand, 0, RI, ...
'RelTol',RelTol_integral,'AbsTol',AbsTol_integral);
% psi = integral(integrand, 0, RI, ...
% 'RelTol',RelTol_integral,'AbsTol',AbsTol_integral);
elseif integral_method == 1
R = linspace(0,RI,nr);
psi = trapz(R,integrand(R));
end
U(2*M+2*N+q) = -(1i*b1/omega(idx)) * (2/(RI^2 * dbesselj(l, chi(q)*RI))) * psi;
end
if solution_method == -3
Unr = U(2*M+N+1:2*M+2*N+Q) - S(2*M+N+1:2*M+2*N+Q,1:2*M+N)*Tr;
Tnr = S(2*M+N+1:2*M+2*N+Q,2*M+N+1:2*M+2*N+Q)\Unr;
fT = Tnr;
elseif solution_method == 3
fT = S(2*M+N+1:2*M+2*N+Q,1:2*M+N)*Tr + ...
S(2*M+N+1:2*M+2*N+Q,2*M+N+1:2*M+2*N+Q)*T - ...
U(2*M+N+1:2*M+2*N+Q);
end
fY = fT;
elseif solution_method == -4 || solution_method == 4
T = Y;
a_vec = T(M+1:M+N);
b_vec = T(1:M);
c_vec = T(M+N+1 : 2*M+N);
d_vec = T(2*M+N+1:2*M+2*N);
e_vec = T(2*M+2*N+1:end);
psi = zeros(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;
if integral_method == 0
psi(q) = quadgk(integrand, 0, RI, ...
'RelTol',RelTol_integral,'AbsTol',AbsTol_integral);
% psi(q) = integral(integrand, 0, RI, ...
% 'RelTol',RelTol_integral,'AbsTol',AbsTol_integral);
elseif integral_method == 1
R = linspace(0,RI,nr);
psi(q) = trapz(R,integrand(R));
end
end
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
if solution_method == -4
fT = S\U;
elseif solution_method == 4
fT = S*T - U;
end
fY = fT;
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
function [Jac,fY,a_vec,b_vec,c_vec,d_vec,e_vec] = JacMatrix(f,Y)
% Forward differencing
[fY,a_vec,b_vec,c_vec,d_vec,e_vec] = f(Y);
dimY = numel(Y);
Jac = zeros(dimY);
for j = 1:dimY
Ysave = Y(j);
% Choose deltaY as to keep increment in the direction of Y(j)
if Ysave == 0
deltaY = sqrt(eps*1e-5);
else
deltaY = sqrt(eps)*sqrt(max(1e-5,abs(Ysave)))*Ysave/abs(Ysave);
end
%% Stepsize control by Hairer/Wanner in ODE integrator RADAU5
%deltaY = sqrt(eps*max(1e-5,abs(Ysave)));
%% Stepsize control in solver hybrd for systems of nonlinear equations
%deltaY = sqrt(eps)*max(abs(Ysave),1)
Y(j) = Y(j) + deltaY;
fYvar = f(Y);
Y(j) = Ysave;
Jac(:,j) = (fYvar - fY)/deltaY;
end
% Central differencing
%
%dimY = numel(Y);
%Jac = zeros(dimY);
%for j = 1:dimY
% Ysave = Y(j);
% if Ysave == 0
% deltaY = sqrt(eps*1e-5);
% else
% deltaY = sqrt(eps)*sqrt(max(1e-5,abs(Ysave)))*Ysave/abs(Ysave);
% endif
% %deltaY = sqrt(eps*max(1e-5,abs(Ysave)));
% %deltaY = sqrt(eps)*max(abs(Ysave),1)
% Y(j) = Y(j) + deltaY;
% fYvarp = compute_fY(Y,S,U,M,N,Q,Z0d,wk,RI,l,Znd,k,chi,b1,omega,idx,...
% solution_method,...
% integral_method,RelTol_integral,AbsTol_integral,nr);
% Y(j) = Ysave;
% Y(j) = Y(j) - deltaY;
% fYvarm = compute_fY(Y,S,U,M,N,Q,Z0d,wk,RI,l,Znd,k,chi,b1,omega,idx,...
% solution_method,...
% integral_method,RelTol_integral,AbsTol_integral,nr);
% Y(j) = Ysave;
% for i = 1:dimY
% Jac(i,j) = (fYvarp(i) - fYvarm(i))/(2*deltaY);
% end
%end
end
EDIT:
At the beginning of the code, change
% Initialize solution vectors
if solver_method == 11 || solver_method == 21
if abs(solution_method) <= 2
solini = zeros(Q,1);
elseif abs(solution_method) == 3
solini = zeros(N+Q,1);
elseif abs(solution_method) == 4
solini = zeros(2*N+2*M+Q,1);
end
elseif solver_method == 12 || solver_method == 22
if abs(solution_method) <= 2
solini = zeros(2*Q,1);
elseif abs(solution_method) == 3
solini = zeros(2*(N+Q),1);
elseif abs(solution_method) == 4
solini = zeros(2*(2*N+2*M+Q),1);
end
end
to
% Initialize solution vectors
if solver_method == 11 || solver_method == 21
if solution_method == 0
solini = zeros(Q,1);
elseif abs(solution_method) == 1 || abs(solution_method) == 2
solini = zeros(Q,1);
elseif abs(solution_method) == 3
solini = zeros(N+Q,1);
elseif abs(solution_method) == 4
solini = zeros(2*N+2*M+Q,1);
end
elseif solver_method == 12 || solver_method == 22
if solution_method == 0
solini = zeros(Q,1);
elseif abs(solution_method) == 1 || abs(solution_method) == 2
solini = zeros(2*Q,1);
elseif abs(solution_method) == 3
solini = zeros(2*(N+Q),1);
elseif abs(solution_method) == 4
solini = zeros(2*(2*N+2*M+Q),1);
end
end
Javeria
am 24 Dez. 2025
First change
options = optimset('Display','none','TolX',TolX_nleq,'TolFun',TolFun_nleq,...
'MaxIter',MaxIter_nleq,'MaxFunEvals',MaxFunEvals_nleq,...
'AutoScaling',AutoScaling_nleq,'ComplexEqn',ComplexEqn_nleq,...
'Updating',Updating_nleq,'TypicalX',TypicalX_nleq,...
'Jacobian',Jacobian_nleq);
to
% options = optimset('Display','none','TolX',TolX_nleq,'TolFun',TolFun_nleq,...
% 'MaxIter',MaxIter_nleq,'MaxFunEvals',MaxFunEvals_nleq,...
% 'AutoScaling',AutoScaling_nleq,'ComplexEqn',ComplexEqn_nleq,...
% 'Updating',Updating_nleq,'TypicalX',TypicalX_nleq,...
% 'Jacobian',Jacobian_nleq);
options = optimset('Display','none','TolX',TolX_nleq,'TolFun',TolFun_nleq,...
'MaxIter',MaxIter_nleq,'MaxFunEvals',MaxFunEvals_nleq);
Now you should also be able to work with MATLAB's "fsolve".
Test this with a small runcase for the following settings:
solution_method = 2;
solver_method = 22;
% Set integration method
integral_method = 0; % 0: Use integral, 1: Use trapezoidal rule
RelTol_integral = 1e-6; % if integral_method = 0: Relative error tolerance in integral solving (default value: 1e-6)
AbsTol_integral = 1e-10; % if integral_method = 0: Absolute error tolerance in integral solving (default value: 1e-10)
nr = 300; % if integral_method = 1: number of discretization points in r-direction
% Set computational 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);
stepX_min = stepX/2^15;
Since it's not certain whether there is one solver which works for all your applications, I don't want to shorten the code. But in essence, there are only 4 settings that you can play with:
solution_method (recommended is = 1 or = 2)
solver_method (recommended is = 12 or - if "fsolve" works with the example above - = 22)
integral_method (= 0 if you want to use the adaptive MATLAB integrator, = 1 if you want to use the trapezoidal rule). If you choose integral_method = 1, you additionally have to set "nr" to the number of evaluation points in the interval [0 RI] for the trapezoidal method.
If it's too confusing for you, you can write a script and call "main" with your settings:
solution_method = ...;
solver_method = ...;
integral_method = ...
nr = ...;
M = ...;
N = ...;
Q = ...;
main(solution_method, solver_method, integral_method, nr, M, N, Q)
Of course, you have to comment out the respective settings in "main" before for that your values are used and not the ones defined in "main".
You should write intermediate results for eta0 to file, not only at the end when eta0 is computed for all values of XkRE (as it is done in the actual code). Even saving all relevant parameters (T,psi,...) for a restart at the value for X_kRE when the solver potentially crahed is advisable.
But be prepared that running cases for big values for M/N/Q (especially Q) will take a long time.
Shouldn't testing
AmA_vec(:,idx) = A*a_vec;
BmB_vec(:,idx) = B*b_vec;
FmC_vec(:,idx) = F*c_vec;
DmD_vec(:,idx) = D*d_vec;
WmE_vec(:,idx) = W*e_vec;
suffice for this purpose ?
Since I don't use MATLAB at the moment, I'd be interested in runtime information for the cases you solve to compare them with octave performance.
Is it possible for you to provide this information together with the values for solution_method, solver_method, integration_method, nr and the case parameters (N,M,Q,RI,RE,h,d,tau,gamma,l) you used ?
Javeria
am 25 Dez. 2025
Javeria
am 25 Dez. 2025
Thank you for the infos.
I also ran the code under octave and got this result for the M = N = Q = 20 case (with no other settings changed except solver_method = 22 and Updating_nleq = 'on' ) within appr. 6 hours. Unfortunately, the "Updating" option is not available under MATLAB - it speeds things up quite significantly, but makes the computation slightly more unstable.
Here is a comparison with the M = N = Q = 10 case:
A1 = readmatrix('data_T_20_20_20_fein_full.txt');
A2 = readmatrix('data_T_10_10_10_fein.txt');
hold on
plot(A1(:,1),A1(:,2))
plot(A2(:,1),A2(:,2))
hold off
grid on
I'm curious and take your case to see how long it takes under octave.
Javeria
am 25 Dez. 2025
At present, for larger problems, I use
solution_method = 2
solver_method = 22
Update_nleq = 'on'
integration_method = 0
with octave.
MATLAB does not have the "Update" option which makes the code significantly faster (but a bit more unstable). And I cannot test MATLAB's "fsolve" and see whether you won't come into trouble with its "exitflag". In the part of the code I wrote, exitflag = 1 indicates that the computation can progress to the next grid point in X_kRE. If exitflag ~= 1, the code tries to insert intermediate points in X_kRE to see whether it can now successfully continue. This can cost quite a lot of time for larger problems and can even make the code stop. And I don't know whether MATLAB's "fsolve" reliably returns exitflag = 1 as the "fsolve" solver under "octave" does (the two solvers are not identical although both share the same name).
Further, until now, I didn't test whether integral_method = 1 makes the code faster. This would also be an option you could test.
So as already adviced above:
Try whether "fsolve" without the Update option (solver_method = 22 instead of solver_method = 12) is still faster than my primitive Newton method code and doesn't make the stepsize being reduced because it doesn't return exitflag = 1. And test integral_method = 1 with an appropriate value for "nr" instead of integral_method = 0.
Or at least test octave instead of MATLAB.
Writing a numerical code takes some time. But most of the time is spent with testing and finding the best solution method.
Torsten
am 25 Dez. 2025
Does your computer and your MATLAB licence admit parallel computing ?
Until now, we use the solution for eta0 in point X_kRE(i-1) as initial guess for the solution in point X_kRE(i) because it stabilizes the computation. So this version of the code would not be parallizable.
But if the solver succeeded also if you could solve for eta(i) independent of what eta0 was in X_kRE(i-1), you could compute in parallel.
Javeria
am 25 Dez. 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.
William Rose
am 31 Dez. 2025
0 Stimmen
10 Kommentare
I don't understand why you think
psi(q) = integral(integrand, 0,r_min, 'AbsTol',1e-8, 'RelTol',1e-6);
is approximately
psi(q) = integral(integrand, 0,RI, 'AbsTol',1e-8, 'RelTol',1e-6);
?
Don't you want to replace
psi(q) = integral(integrand, 0,RI, 'AbsTol',1e-8, 'RelTol',1e-6);
by
psi(q) = integral(integrand, r_min ,RI, 'AbsTol',1e-8, 'RelTol',1e-6);
here ?
Further, you can avoid "integral" calling your function "debug_integrand" with r being a vector with more than one element by using 'Arrayvalued',true as further option when calling "integral":
value_integral = integral(fun,a,b,'AbsTol',...,'RelTol',...,'ArrayValued',true);
This would turn your recursive calls of v_D in this part
% --- Allow vector r during integration ---
if numel(r) > 1
v_D_val = zeros(size(r));
for i = 1:numel(r)
v_D_val(i) = v_D(N, Q, r(i), Z0d, wk, RI, l, d_vec, Znd, k, e_vec, chi);
end
return
end
unnecessary.
Did you validate your results with the results of one of the Newton solvers ?
Here is a plot of your method against the Newton solver in the critical region:
A = load('Results_Newton.txt');
B = load('Results_Iteration.txt');
hold on
plot(A(:,1),A(:,2))
plot(B(:,1),B(:,2))
hold off
grid on
Did you check how large M, N and Q have to be for the parameter set you use at the moment such that your series can be considered as converged ?
Javeria
am 5 Jan. 2026
I used integral_method = 1, nr = 160 , M = N = Q = 60 and the MATLAB version of the solver "nleq1" given at the below link and it took about 11 hours to run under "octave" for 160 output points. You can multiply the 11 hours by approximately 2/3 for MATLAB, I guess, so that you arrive at approximately 7 hours. The solution for eta0 stabilized at these values.
Also using "fsolve" gave acceptable run times of about 15 hours under "octave", thus approximately 10 hours with MATLAB.
Javeria
am 5 Jan. 2026
I think the settings in the attached code are best-possible for your case.
You will have to download the nleq1 code from the website because I'm not allowed to distribute it. Only copy the nleq1.m file - the other files are example files that are unnecessary.
Javeria
am 5 Jan. 2026
Torsten
am 5 Jan. 2026
Are you angry or what do your upper case letters mean ?
Javeria
am 5 Jan. 2026
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.
6 Kommentare
We are still in search of a reliable and fast solution for the case where fixed point iteration cannot work. Do you have a suggestion ?
If it had been obvious right from the beginning that the project would become computationally intensive, I had chosen C or Fortran. Most probably only 1/10 of the computational time needed compared to MATLAB.
Javeria
am 10 Jan. 2026
The software codes provided to solve the nonlinear part of your system of equations (my own simple Newton method, "fsolve" from MATLAB and "nleq1" from ZIB) are general-purpose codes. That means that they are not adopted in any way to solve your special problem for the Q modes. And for this kind of codes, there is nothing else but the formulation of the system of equations itself (thus the f in f(x) = 0) that you can use to influence their performance.
I don't know exactly what you mean by "estimation of damping" as input. I just took your equations and handed them to the solver without backgound knowledge about the physics of your problem. But if you intend to use "estimation of damping" as a means to influence the treatment of the nonlinearity and/or to accelerate convergence, my guess is that it won't be possible in the framework of these general purpose codes. To take advantage of special structures or features of your system to be solved, I guess you will have to implement it in your own nonlinear solver.
Javeria
am 12 Jan. 2026
Javeria
am 12 Jan. 2026
32 Kommentare
The way you compute the integral is important for the runtime of your code, but is not responsible for and can't avoid the divergence. The fixed-point iteration for T in the part
psi = zeros(Q,1);
for iter = 1:max_iter
% (A) Update 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
% (B) Solve the linear system
T = S \ U;
% (C) Extract coefficients
b_vec = T(1:M);
a_vec = T(M+1 : M+N);
c_vec = T(M+N+1 : 2*M+N);
d_vec = T(2*M+N+1:2*M+2*N);
e_vec = T(2*M+2*N+1:end);
% Somehow solve the integrals here
psi = something
% (E) Check convergence
if ~isempty(T_old)
diff_T = max(abs(T - T_old));
if diff_T < tol
converged = true;
break
end
end
T_old = T;
end % End iteration loop
diverges because the fixed-points for your second parameter case are non-attractive. This results in unrealistic values for d_vec = T(2*M+N+1:2*M+2*N) and e_vec = T(2*M+2*N+1:end) in the integration with the consequence that the integration also returns Inf or NaN values. Thus it won't help to modify the integration - it is ok with "integral" or "trapz" as is.
You will have to work on the method how to solve
S*T = U(T)
or
S*T-U(T) = 0
At the moment, you iterate T as
T_(n+1) = S\U(T_n) (1)
As I tried to explain in one of my comments from above (and as I proved with the code I use), the function
g(T) = S\U(T)
has eigenvalues of absolute value > 1 in the solution points T* of your system S*T - U(T) = 0 for certain parameter constellations. Thus when you use the fixed-point iteration for these parameters, you will never get T*.
Just to show you what I mean: For your code example from above, I get the following maximum absolute values of the eigenvalues of the iteration function in the solution point T* for the 10 values of X_kRE:
4.4e-3, 0.895, 6.8785, 19.356, 5.3709, 2.9084, 1.7406, 1.1185, 0.7607, 0.5418.
As you can see, from theory, you should only get convergence of the fixed-point method for the first and last two X_kRE values. You were lucky that the third last iteration converged, too.
Sometimes modifying the iteration function can solve such problems. E.g. if you want to solve
tan(x) = x
and iterate
x_(n+1) = tan(x_n)
you get divergence, but if you set x = atan(y) and write your equation as
y = atan(y)
and iterate
y_(n+1) = atan(y_n)
your fixed-point iteration converges (because the iteration function g(y) = atan(y) has derivative g'(y) = 1/(1+y^2) which is < 1).
Maybe you can find such a transformation for your case. Do you have an idea how exactly U(2*N+2*M+1:2*N+2*M+Q) look like as functions of T ? Are they quasi-quadratic in T ?
I didn't think about it more deeply and simply applied Newton's method to compute e_vec = T(2*M+2*N+1:end). For the other components of T you have linear equations and thus a_vec, b_vec, c_vec and d_vec can be deduced from e_vec.
I will have a look in the article later.
Torsten
am 12 Jan. 2026
I read the article and - correct me if I'm wrong - the problem considered in the article is instationary while yours is stationary. Thus the numerical tool used in the article is a solver for ordinary differential equations while your tool must be a solver for algebraic equations. I don't see how to extract useful information for your case.
Javeria
am 13 Jan. 2026
the only difference is that they have time dependent study and i have the time independent they also used the same method and form the system of equations and for the non linear integral they mentioned that method.
I know this. The authors of the article will most probably use a scientific code that solves their system of ordinary differential equations. In MATLAB, this would be "ode15s". For your problem, MATLAB offers "fsolve" to solve nonlinear systems of algebraic equations. It's unusual that authors develop their own solution method (like you try to do) - they use general-purpose software because it's mostly faster and more reliable. And concerning the integral: you can use the author's method if you like, but it won't change anything in the behaviour of your fixed-point iteration.
Let me put it clear:
Newton's method is the standard method to solve systems of nonlinear equations. If you find an article that solves some kind of stationary problem like you do, you will find with 99% probability that the authors used Newton's method. And it works as you can see from the code I use. If you only want to use it for your own work and don't want to sell it, speed doesn't matter in my opinion as long as you get results in an acceptable time. For small values for Q, you get results almost instantly. If you run it for larger values, you can start several jobs in parallel and get the results next day. For me, this is acceptable. If you further want to speed up the code, start anew and use Fortran or C.
If you have questions concerning the code or an idea how to use another method instead of your fixed-point iteration, you are again welcome to ask.
Javeria
am 13 Jan. 2026
The code I provided uses Newton-Raphson to solve for e_vec and Gauss-Kronos quadrature to determine the integral. Why do you think Newton-Raphson for (24) (whatever it is) and 6-point Gaussian quadrature will be better ?
If you find references that use Newton-Raphson, the respective code will behave similarly in speed as the one I provided.
Paul
am 13 Jan. 2026
"the function
g(T) = S\U(T)
has eigenvalues of absolute value > 1"
If the iteration function is
T_(n+1) = S\U(T_n) (1)
then shouldn't we be concerned about eigenvlaues of
inv(S)
rather than
g(T) = S\U(T) ?
I'm sure that's how it should be if U(T) were simply U(T) = T.
Sorry, I meant the eigenvalues of the Jacobian of the iteration function
g(T) = S\U(T)
with respect to T here.
If U(T) = T, this will be
inv(S)
If U is a general function of T, it should be
inv(S) * Jac(U)
I guess.
Javeria
am 13 Jan. 2026
Newton-Raphson is used for solution_method = 1, 2, 3 or 4.
Gauss-Kronos quadrature is used for integral_method = 1.
Solver_method = 12, 22 or 32 can be selected.
I recommend solution_method = 2, integral_method = 2 (the trapoidal rule gave almost identical results to Gauss-Kronos, but was a lot faster for bigger problems) and solver_method = 32 (nleq1 was the fastest of the Newton-solvers).
I don't want to destroy choices because I'm a pessimist and want to have a plan B if plan A doesn't succeed.
The abscissa for plotting is computed as RE * X_kRE * tanh(X_kRE * h). I saw that you changed this to
2*pi / sqrt(g * X_kRE * tanh(X_kRE * h)) in your last posted code.
Javeria
am 14 Jan. 2026
The blue curve is the curve that is created at the X_kRE values you prescribe.
The solution is advanced from one X_kRE value to the next, and the initial guess for the vector of unknowns at step i is taken as the solution at step (i-1). Sometimes it may happen that the solver doesn't converge at step i. Then a value for your X_kRE is inserted in the middle of the interval [X_kRE(i-1),X_kRE(i)] and the solver is called for (X_kRE(i-1)+X_kRE(i))/2, again with the solution at step (i-1) as initial guess. The red curve shows the solution in your prescribed values for X_kRE plus the solution in the points for X_kRE that were inserted.
I write two files to your working directory: data_test and data_test_full. data_test contains the solution in your prescribed X_kRE vector, data_test_full contains the solution in your prescribed X_kRE vector plus in the points that were inserted. You can use these files later to make plots for different cases. Don't forget to rename the files after a run has finished - otherwise they will be overwritten in the next run.
Further there is an option to restart the computation from the results obtained at a certain value for X_kRE if a run takes too long and you want to shut down your computer. The name of the restart file is solution_test.
You can prescribe the names of all the files written in the section "Solution file names" of the code.
I understand that you want to stick to your own coding. I also thought when I started working in the field of numerical analysis: Why do people need so much lines of code for this problem ? I can do it equally well with 20 lines. I was wrong.
Your Newton method is theoretically fine. I also started to use the complex Newton iteration as you do, but soon asked myself: Why does the solver need so many iterations ? It turned out that the Jacobian matrix was dependent on the direction in which I computed the partial derivatives. So if you choose eps_fd = 1e-6*i (i = complex unit), e.g., you will get another J. This shouldn't happen if you want to apply complex Newton iteration successfully - the results for J shouldn't depend on the direction of eps_fd. What is the reason ? The function you integrate is not holomorphic. As soon as I started to set up the problem as real(f) = 0 and imag(f) = 0 - thus the usual real Newton iteration applied to real and imaginary part of f - this worked much better. And if you increase the value of Q, you will soon see that the "primitive" Newton methods you and I coded are much too slow.
But the main reason for the discrepanies between your and my results (in my opinion - I didn't test it) (see below) is that you only use 6 points to discretize the integral [0 RI]. Although RI is small here, this cannot be a sufficient resolution. 6-point quadrature doesn't mean that - independent of the length of the interval - you only use 6 points in which you evaluate f. It means to subdivide the complete interval into a certain number of subintervals, apply the stencil to each of these subintervals and sum the results.
So for the time to come, my advice is to at least compare results from your solution with the solutions from the solver I provided before proceeding.
A1 = dlmread('Results_Javeria.txt');
A2 = dlmread('Results_Torsten.txt');
hold on
plot(A1(:,1),A1(:,2),'r');
plot(A2(:,1),A2(:,2),'b');
grid on
hold off
xlim([1 1.75])
Javeria
am 15 Jan. 2026
I think increasing N will not make a big difference - increasing Q is the key factor for runtime. But I didn't test increasing N alone while keeping M and Q constant.
This code is almost equivalent to the settings
solution_method = 1
solver_method = 11
integral_method = 2
nr = 160
in the code I use and gives similar results:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% NEWTON–RAPHSON NONLINEAR SOLVER
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% --- 6-point Gauss–Legendre quadrature on [0, RI]
%xg = [-0.9324695142; -0.6612093865; -0.2386191861; ...
% 0.2386191861; 0.6612093865; 0.9324695142];
%wg = [ 0.1713244924; 0.3607615730; 0.4679139346; ...
% 0.4679139346; 0.3607615730; 0.1713244924];
%
%r_g = 0.5 * RI * (xg + 1);
%w_g = 0.5 * RI * wg;
% --- Newton parameters
psi = zeros(Q,1); % initial guess (linear solution)
tol_NR = 1e-6;
max_NR = 100;
nr = 160;
R = linspace(0,RI,nr);
for it = 1:max_NR
% ===============================
% 1) Evaluate F(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
% Solve linear system
T = S \ U;
d_vec = T(2*M+N+1 : 2*M+2*N);
e_vec = T(2*M+2*N+1 : end);
% Compute nonlinear residual F
Fnl = zeros(Q,1);
for q = 1:Q
%Iq = 0;
%for kq = 1:6
% r = r_g(kq);
% v = v_D(N,Q,r,Z0d,wk,RI,l,d_vec,Znd,k,e_vec,chi);
% Iq = Iq + w_g(kq) * abs(v) * v ...
% * besselj(l,chi(q)*r) * r;
%end
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;
Fnl(q) = psi(q) - trapz(R,integrand(R));
end
% Convergence check
if norm(Fnl,inf) < tol_NR
fprintf('Newton converged in %d iterations\n', it);
break
end
% ===============================
% 2) Jacobian (finite difference)
% ===============================
J = zeros(Q,Q);
eps_fd = 1e-6;
for j = 1:Q
psi_p = psi;
psi_p(j) = psi_p(j) + eps_fd;
for q = 1:Q
U(2*M+2*N+q) = ...
-(1i*b1/omega(idx)) ...
* (2/(RI^2*dbesselj(l,chi(q)*RI))) * psi_p(q);
end
T_p = S \ U;
d_p = T_p(2*M+N+1 : 2*M+2*N);
e_p = T_p(2*M+2*N+1 : end);
for q = 1:Q
%Iqp = 0;
%for kq = 1:6
% r = r_g(kq);
% v = v_D(N,Q,r,Z0d,wk,RI,l,d_p,Znd,k,e_p,chi);
% Iqp = Iqp + w_g(kq) * abs(v) * v ...
% * besselj(l,chi(q)*r) * r;
%end
integrand = @(r) abs(v_D(N,Q,r,Z0d,wk,RI,l,d_p,Znd,k,e_p,chi)) ...
.* v_D(N,Q,r,Z0d,wk,RI,l,d_p,Znd,k,e_p,chi) ...
.* besselj(l, chi(q)*r) .* r;
J(q,j) = (psi_p(q) - trapz(R,integrand(R)) - Fnl(q)) / eps_fd;
end
end
% ===============================
% 3) Newton update
% ===============================
delta = J \ Fnl;
psi = psi - delta;
end
if it == max_NR
warning('Newton did not converge');
end
Javeria
am 24 Jan. 2026
For the first question:
we need your code to reproduce the problem.
For the second question:
b1 only appears in U(2*M+2*N+q) as a multiplicative factor. You could try to compute the sensitivity of eta0 with respect to b1 (i.e. d(eta0)/d(b1)). Maybe in your model, its influence on eta0 is neglectable small.
Changing b1 gives a different height of the peak for eta0. Further, the curve for eta0 is slightly translated to the left.
A = dlmread('data_b1=normal.txt');
AA = dlmread('data_b1_with_H=100.txt');
hold on
plot(A(:,1),A(:,2),'r')
plot(AA(:,1),AA(:,2),'b')
hold off
grid on
Javeria
am 25 Jan. 2026
I can't help you with the question whether your model is physically adequate or not.
I can only look for mathematical problems.
Your integration with the 6-point gauss-Legendre quadrature is still not corrected and a corrected code gives different peak heights for eta0 if b1 is varied.
Javeria
am 26 Jan. 2026
I corrected your code with "trapz" instead of your 6-point Gauss-Legendre quadrature in this discussion. Alternatively, you can use MATLAB's "integral" as before. If you sufficiently resolve the region around the peak for "eta0" (I used 200 instead of 20 points), you should see the different heights for different values of H also with your code. If it doesn't converge for certain values of X_kRE, try my version (I got the results with my version).
Javeria
am 27 Jan. 2026
did you meant just replace my 6-point Gauss-Legendre quadrature by "trapz" as following
Yes, or by
Fnl(q) = psi(q) - integral(integrand, 0, RI, 'RelTol',1e-6,'AbsTol',1e-10);
resp.
J(q,j) = (psi_p(q) - integral(integrand, 0, RI, 'RelTol',1e-6,'AbsTol',1e-10) - Fnl(q)) / eps_fd;
I'd use the version with "integral" because the choice of the resolution of [0 RI] = [0 34.5] with nr = 160 points is quite arbitrary and it's not obvious that it is sufficient for a good approximation of the integral.
"integral" uses an internal estimate of the error and adjusts the number of points as to match your given error tolerances RelTol and AbsTol.
Javeria
am 28 Jan. 2026 um 2:22
Javeria
am 28 Jan. 2026 um 4:12
If the integral is computed reliably, it can make it harder for Newton's method to converge. I saw this especially in the peak region for eta0. But this doesn't mean that you should change the way you compute the integral, but the way you coded Newton's method. Using MATLAB's "integral" is the best you can do (although it's slower than "trapez") because it has an internal error control (thus it uses as many quadrature points as needed to satisfy your given tolerable error constants RelTol and AbsTol).
Here are some tricks and necessary changes I used for Newton's method in the code I provided:
Often, a constant initial guess for psi (psi = zeros(Q,1)) in all steps is not advisable.That's why in the code I provided, I use the psi I get in iteration (i-1) as start point in iteration i.
If this still does not make Newton's method converge, I half the distance between X_kRE(i-1) and X_kRE(i) and try to compute psi at (X_kRE(i-1) + X_kRE(i))/2 - again with psi(i-1) as initial guess.
And I rewrote the problem as real(f) = 0 and imag(f) = 0 because complex Newton iteration for non-holomorphic f causes problems in the computation of the Jacobian.
Javeria
am 28 Jan. 2026 um 12:17
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
am 29 Jan. 2026 um 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).
Kategorien
Mehr zu Linear Algebra finden Sie in Hilfe-Center und File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!





















