numerical Instabilities for bessel functions

322 visualizzazioni (ultimi 30 giorni)
Javeria
Javeria il 8 Dic 2025
Modificato: Torsten il 15 Gen 2026 alle 21:35
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 Commenti
Javeria
Javeria il 8 Dic 2025
@Walter Roberson@Sam Chak what i have that we need to treat this non linear term iteratively to solve our system of equation;
The principle of iteration is the same. You may go in another way: 1) assuming your coefficient b1=0, you should have everything including v_D(r);
2) using a new a’( r)=a+b|v_D( r)| in place of a to solve all equations to get a new v_D( r); 3) to check the new v_D with previous one, if the difference is larger than your tolerance, continue the step 2).
David Goodmanson
David Goodmanson il 23 Dic 2025 alle 10:25
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.

Accedi per commentare.

Risposte (4)

Torsten
Torsten il 9 Dic 2025
Modificato: Torsten il 9 Dic 2025
Although it's slow, the full Newton method seems to work in this case (except for one value of X_kRE):
clear all;
close all;
clc;
tic;
% diary('iter_log.txt');
% diary on
% fprintf('Run started: %s\n', datestr(now));
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Parameters from your setup
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
old = false;
if old
N = 20; % Number of N modes
M = 20; % Number of M modes
Q = 20; % Number of Q modes
RI = 34.5; % Inner radius (m)
RE = 47.5; % Outer radius (m)
h = 200; % Water depth (m) % Interface depth (m)
d = 38; % Interface depth (m) % Draft (m)
tau = 0.1; % porosity ratio
X_kRE = linspace(0.0025, 0.16,100); %%%%%% this is now k_0
else
N = 10; % Number of N modes
M = 10; % Number of M modes
Q = 10; % Number of Q modes
RI = 0.4; % Inner radius (m)
RE = 0.8; % Outer radius (m)
h = 1.0; % Water depth (m) % Interface depth (m)
d = 0.3; % Interface depth (m) % Draft (m)
tau = 0.2; % porosity ratio
X_kRE = linspace(0.05, 4.5, 40); %%%%%% this is now k_0
end
b = h-d; % Distance from cylinder bottom to seabed (m)
g = 9.81; % Gravity (m/s^2)
gamma = 1.0; % discharge coefficient
b1 = (1 - tau) / (2 * gamma * tau^2); % nonlinear coeff
tol = 1e-4; % convergence tolerance
max_iter = 100; % max iterations
l = 0; % Azimuthal order
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% === Compute Bessel roots and scaled chi values ===
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
roots = bessel0j(l, Q); % Zeros of J_l(x)
chi = roots ./ RI; % chi_q^l = roots scaled by radius
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% ==============Initialize eta0 array before loop
%===============Initialize omega array before loop
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
eta0 = zeros(size(X_kRE)); % Preallocate
omega = zeros(size(X_kRE)); %%%% omega preallocate
% iters_used = zeros(length(X_kRE),1);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for idx = 1:length(X_kRE)
wk = X_kRE(idx); % dimensional wavenumber
omega(idx) = sqrt(g * wk * tanh(wk * h)); % dispersion relation
% fprintf('\n--- idx %3d (T=%6.3f s) ---\n', idx, 2*pi/omega(idx));
% drawnow;
%%%%%%%%%%%%% Compute group velocity %%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Cg(idx) = (g * tanh(wk * h) + g * wk * h * (sech(wk * h))^2) ...
* omega(idx) / (2 * g * wk * tanh(wk * h));
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% =======Compute a1 based on tau
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
a1(idx) = 0.93 * (1 - tau) * Cg(idx);
% dissip(idx) = a1(idx)/omega(idx);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%=============== Find derivative of Z0 at z = -d
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Z0d = (cosh(wk * h) * wk * sinh(wk * b)) / (2 * wk * h + sinh(2 * wk * h));
fun_alpha = @(x) omega(idx)^2 + g*x.*tan(x*h); %dispersion equation for evanescent modes
opts_lsq = optimset('Display','none'); % ← add this
%opts_lsq = [];
for n = 1:N
if n == 1
k(n) = -1i * wk;
else
x0_left = (2*n - 3) * pi / (2*h);
x0_right = (2*n - 1) * pi / (2*h);
guess = (x0_left + x0_right)/2;
% Use lsqnonlin for better convergence
k(n) = lsqnonlin(fun_alpha, guess, x0_left, x0_right,opts_lsq);
%%%%%%%%%%%%% Derivative of Z_n at z = -d %%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Znd(n) = -k(n) * (cos(k(n)*h) * sin(k(n)*b)) / (2*k(n)*h + sin(2*k(n)*h));
end
end
%% finding the root \lemda_m =m*pi/b%%%%%%
for j=1:M
m(j)=(pi*(j-1))/b;
end
% %%%%%%%%%%%%%%%%%%%%% Define matrix A %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
A = zeros(M,N);
A(1,1) = (cosh(wk*h)*sinh(wk*b))/(wk*b*(2*wk*h+sinh(2*wk*h)))*(besselh(l, wk*RE)) /(dbesselh(l, wk*RE)); % Set A_00
for j = 2:N
A(1,j) =(cos(k(j)*h)*sin(k(j)*b) )/ (k(j)*b*(2*k(j)*h + sin(2*k(j)*h)))*( besselk(l, k(j)*RE) )...
/(dbesselk(l, k(j)*RE)); % Set A_0n
end
for i = 2:M
A(i,1) =(2 * cosh(wk*h) * (-1)^(i-1) * wk*sinh(wk*b)) / (b * (2*wk*h + sinh(2*wk*h))*(wk^2 + m(i)^2))...
* (besselh(l, wk*RE))/ (dbesselh(l, wk*RE)); % Set A_m0
end
for i = 2:M
for j = 2:N
A(i,j) = (2*cos(k(j)*h)*(-1)^(i-1) *k(j)* sin(k(j)*b))/ (b * (2*k(j)*h + sin(2*k(j)*h))*(k(j)^2-m(i)^2))...
*(besselk(l, k(j)*RE)) / (dbesselk(l, k(j)*RE));% Set A_mn
end
end
%%%%%%%%%%%%%%%%%%%%%%%% DefineMatrix B %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
B=zeros(N,M);
B(1,1) = (4*sinh(wk*b)) / (RE * wk*log(RE/RI) *cosh(wk*h)); %set B_00
%%%%%%%%%%%B_0m%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for j = 2:M
Rml_prime_RE = m(j)*(besselk(l, m(j)*RI) * dbesseli(l, m(j)*RE) ...
- besseli(l, m(j)*RI) * dbesselk(l, m(j)*RE))...
/(besselk(l, m(j)*RI) * besseli(l, m(j)*RE) - besselk(l, m(j)*RE) * besseli(l, m(j)*RI));
B(1,j) = Rml_prime_RE * (4 * wk * (-1)^(j-1) * sinh(wk*b)) / (cosh(wk*h) * (wk^2 + m(j)^2));
end
%%%%%%%%%%%B_n0%%%%%%%%%%%%%%%%%%%%%%%%%%%
for i=2:N
B(i,1)=(4*sin(k(i)*b))/(RE*k(i)*log(RE/RI)*cos(k(i)*h));% Set B_n0
end
%%%%%%%%%%%%%%%%%%%%%%%%%%B_nm%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for i=2:N
for j=2:M
Rml_prime_RE = m(j)*(besselk(l, m(j)*RI) * dbesseli(l, m(j)*RE) ...
- besseli(l, m(j)*RI) * dbesselk(l, m(j)*RE))...
/(besselk(l, m(j)*RI) * besseli(l, m(j)*RE) - besselk(l, m(j)*RE) * besseli(l, m(j)*RI));
B(i,j) = Rml_prime_RE * (4 * k(i) * (-1)^(j-1) * sin(k(i)*b)) / (cos(k(i)*h) * (k(i)^2 - m(j)^2));
end
end
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Define Matrix C %%%%%%%%%%%%%%%%%%%%%%%%
C = zeros(N,M);
C(1,1) = -4 * sinh(wk*b) / (RE * wk*log(RE/RI) * cosh(wk*h));
% %%%%% DEfine C0m%%%%%%
for j = 2:M
Rml_prime_star_RE = m(j)*(besseli(l, m(j)*RE) * dbesselk(l, m(j)*RE) ...
- besselk(l, m(j)*RE) * dbesseli(l, m(j)*RE))...
/(besselk(l, m(j)*RI) * besseli(l, m(j)*RE) - besselk(l, m(j)*RE) * besseli(l, m(j)*RI));
C(1,j) = Rml_prime_star_RE * (4 * wk * (-1)^(j-1) * sinh(wk*b)) / (cosh(wk*h) * (wk^2 + m(j)^2));
end
% %%%%%%% Define Cn0%%%%%%
for i = 2:N
C(i,1) = -4 * sin(k(i)*b) / (RE *k(i)* log(RE/RI) * cos(k(i)*h));
end
% % %%%%%% Define Cnm%%%%%%
for i = 2:N
for j =2:M
Rml_prime_star_RE = m(j)*(besseli(l, m(j)*RE) * dbesselk(l, m(j)*RE) ...
- besselk(l, m(j)*RE) * dbesseli(l, m(j)*RE))...
/(besselk(l, m(j)*RI) * besseli(l, m(j)*RE) - besselk(l, m(j)*RE) * besseli(l, m(j)*RI));
C(i,j) = Rml_prime_star_RE * (4 * k(i) * (-1)^(j-1) * sin(k(i)*b)) / (cos(k(i)*h) * (k(i)^2 - m(j)^2));
end
end
%%%%%%% write Matrix D %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
D = zeros(M,N);
%%% write D00 %%%%%%
D(1,1) = (cosh(wk*h) * sinh(wk*b)) / (wk*b * (2*wk*h + sinh(2*wk*h))) * (besselj(l, wk*RI))/ (dbesselj(l, wk*RI));
%%%%% write D0n %%%%%%%%%%
for j =2:N
D(1,j) = (cos(k(j)*h) * sin(k(j)*b)) / (k(j)*b * (2*k(j)*h + sin(2*k(j)*h))) * (besseli(l, k(j)*RI) )/ (dbesseli(l, k(j)*RI));
end
%%%%%% write Dm0 %%%%%%%%%
for i = 2:M
D(i,1) = (2 * cosh(wk*h) * (-1)^(i-1) * wk * sinh(wk*b)) /(b * (2*wk*h + sinh(2*wk*h)) * (wk^2 + m(i)^2))...
*(besselj(l, wk*RI) )/(dbesselj(l, wk*RI));
end
%%%%% Define Dmn%%%%%%%%
for i = 2:M
for j = 2:N
D(i,j) = (2 * cos(k(j)*h) * (-1)^(i-1) * k(j) * sin(k(j)*b)) /(b * (2*k(j)*h + sin(2*k(j)*h)) * (k(j)^2 - m(i)^2))...
*(besseli(l, k(j)*RI)) / (dbesseli(l, k(j)*RI));
end
end
%%%%%%%%%%%%%%%%%%%%%%%%%% Define Matrix E %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
E = zeros(N, M); % Preallocate the matrix E
%%%%% Define E00%%%%%%%%%%%%
E(1,1) = (4 * sinh(wk*b)) / (RI *wk* log(RE/RI) * cosh(wk*h));
%%%% Define Eom %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for j = 2:M
Rml_prime_RI = m(j)*(besselk(l, m(j)*RI) * dbesseli(l, m(j)*RI) ...
- besseli(l, m(j)*RI) * dbesselk(l, m(j)*RI))...
/ (besselk(l, m(j)*RI) * besseli(l, m(j)*RE) - besselk(l, m(j)*RE) * besseli(l, m(j)*RI));
E(1,j) = Rml_prime_RI * (4 * wk * (-1)^(j-1) * sinh(wk*b)) / (cosh(wk*h) * (wk^2 + m(j)^2));
%
end
%%%%%% Define Eno%%%%%%%%%%%%%%
for i = 2:N
E(i,1) = (4 * sin(k(i)*b)) / (RI *k(i)* log(RE/RI) * cos(k(i)*h));
end
for i = 2:N
for j = 2:M
Rml_prime_RI = m(j)*(besselk(l, m(j)*RI) * dbesseli(l, m(j)*RI) ...
- besseli(l, m(j)*RI) * dbesselk(l, m(j)*RI))...
/ (besselk(l, m(j)*RI) * besseli(l, m(j)*RE) - besselk(l, m(j)*RE) * besseli(l, m(j)*RI));
E(i,j) = Rml_prime_RI * (4 * k(i) * (-1)^(j-1) * sin(k(i)*b)) / (cos(k(i)*h) * (k(i)^2 - m(j)^2));
end
end
%%%%%%%%%%%%%%%%%%%%%%%% Now write matrix F %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
F = zeros(N,M);
%%%%%%%% Define F00 %%%%%%%%%%
F(1,1) = (-4 * sinh(wk*b)) / (RI * wk*log(RE/RI) * cosh(wk*h));
% %%%%%% Define F0m%%%%
for j = 2:M
Rml_star_prime_RI = m(j)*(besseli(l, m(j)*RE) * dbesselk(l, m(j)*RI) ...
- besselk(l, m(j)*RE) * dbesseli(l, m(j)*RI))/((besselk(l, m(j)*RI) ...
* besseli(l, m(j)*RE) - besselk(l, m(j)*RE) * besseli(l, m(j)*RI)));
F(1,j) = Rml_star_prime_RI * (4 * wk * (-1)^(j-1) * sinh(wk*b)) / (cosh(wk*h) * (wk^2 + m(j)^2));
end
%%%%% Defin Fn0%%%%%%
for i = 2:N
F(i,1) = (-4 * sin(k(i)*b)) / (RI *k(i)* log(RE/RI) * cos(k(i)*h));
end
%%%%%%% Define Fnm %%%%%%%
for i = 2:N
for j = 2:M
Rml_star_prime_RI = m(j)*(besseli(l, m(j)*RE) * dbesselk(l, m(j)*RI) ...
- besselk(l, m(j)*RE) * dbesseli(l, m(j)*RI))/((besselk(l, m(j)*RI) ...
* besseli(l, m(j)*RE) - besselk(l, m(j)*RE) * besseli(l, m(j)*RI)));
F(i,j) = Rml_star_prime_RI * (4 * k(i)* (-1)^(j-1) * sin(k(i)*b)) / (cos(k(i)*h) * (k(i)^2 - m(j)^2));
end
end
% % %%%%%%%%%% Write Matrix W %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
W = zeros(N,Q);
% %%W_{0q}
for q = 1:Q
r = exp(-2*chi(q)*b); % e^{-2*chi*b}
coth_chib = (1 + r) / (1 - r); % coth(chi*b) via decaying exp
W(1,q) = -(4*chi(q)/cosh(wk*h))*(( wk*sinh(wk*b)*coth_chib - chi(q)*cosh(wk*b) ) / ( wk^2 - chi(q)^2 ) ); %prof.chen
%
% % W(1, q) = - (4 * chi(q) / cosh(wk * h)) * (chi(q) * sinh(chi(q) * b) * cosh(wk * b) - wk * sinh(wk * b) * cosh(chi(q) * b)) ...
% % / ((chi(q)^2 - wk^2) * sinh(chi(q) * b));
end
% % %%% Write W_{nq}
for n = 2:N
for q = 1:Q
r = exp(-2*chi(q)*b); % e^{-2*chi*b}
coth_chib = (1 + r) / (1 - r); % coth(chi*b), stable
W(n,q) = -(4*chi(q)/cos(k(n)*h)) * ...
( ( k(n)*sin(k(n)*b)*coth_chib + chi(q)*cos(k(n)*b) ) ... %prof.chen
/ ( k(n)^2 + chi(q)^2 ) );
% % W(n, q) = - (4 * chi(q) / cos(k(n) * h)) * (chi(q) * sinh(chi(q) * b) * cos(k(n) * b) + k(n) * sin(k(n) * b) * cosh(chi(q) * b)) ...
% % / ((chi(q)^2 + k(n)^2) * sinh(chi(q) * b));
end
end
%%%%%%%%%%%%%%%%%%%%%%%Find E1_{q}%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
f2 = omega(idx)^2 / g;
E1 = zeros(Q,1);
S_q = zeros(Q,1);
for q = 1:Q
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% this expression is same as prof. Chen code
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
exp2chi_d = exp(-2 * chi(q) * d); % e^{-2*chi*d}
exp2chi_b = exp(-2 * chi(q) * b); % e^{-2*chi*(h-d)}
exp2chi_h = exp(-2 * chi(q) * h); % e^{-2*chi*h}
exp1chi_d = exp(-chi(q) * d); % e^{-chi*d}
% intermediate ratio (same order as Fortran)
BB = (f2 - chi(q)) * (exp2chi_b / exp2chi_d);
BB = BB + (f2 + chi(q)) * (1 + 2 * exp2chi_b);
BB = BB / (f2 - chi(q) + (f2 + chi(q)) * exp2chi_d);
BB = BB / (1 + exp2chi_b);
% final coefficient
E1(q) = 2 * (BB * exp2chi_d - 1);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% This expression is that one derived directly
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% % N0, D0, H0 exactly as in the math above
% N0 = (chi(q) - f2) + (chi(q) + f2) * exp2chi_h;
% D0 = (f2 - chi(q)) + (f2 + chi(q)) * exp2chi_d;
% H0 = 1 + exp2chi_b;
%
% % final coefficient (no add/subtract step)
% E1(q) = 2 * N0 / (D0 * H0);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% This one is the hyperbolic one form
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% num = chi(q)*cosh(chi(q)*h) - f2*sinh(chi(q)*h);
% den = (f2*cosh(chi(q)*d) - chi(q)*sinh(chi(q)*d)) * cosh(chi(q)*(h - d)); % or: * cosh(chi(q)*b)
% E1(q) = num / den;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%% This is the upper region velocity potential at z=0, write
%%%%%%%%%%%%%%%% from prof.chen code
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
CCmc = 1 + (exp2chi_b / exp2chi_d) + 2*exp2chi_b;
CCmc = CCmc / (1 + exp2chi_b);
CCmc = BB * (1 + exp2chi_d) - CCmc;
S_q(q) = CCmc * exp1chi_d; % this equals the summand value for this q
% --- direct hyperbolic (exact) ---
% % % C_hyp = ( chi(q)*cosh(chi(q)*h) - f2*sinh(chi(q)*h) ) ...
% % % / ( (f2*cosh(chi(q)*d) - chi(q)*sinh(chi(q)*d)) * cosh(chi(q)*b) );
% % %
% % % S_hyp(q) = C_hyp * cosh(chi(q)*d) + sinh(chi(q)*h) / cosh(chi(q)*b); % summand
end
% H(q): diagonal term (exponential form)
H1 = zeros(Q,1);
for q = 1:Q
% dissip =
r = exp(-2*chi(q)*b); % e^{-2*chi*b}
mid = 4*r/(1 - r^2) - E1(q); % 2/sinh(2*chi*b) - C_m
H1(q) = mid + 1i*a1(idx)*chi(q)/omega(idx); % mid + i*(a1/w)*chi
end
H = diag(H1);
%%%%%%%%%%matric G%%%%%%%%%%%%%%%%%%
G = zeros(Q, N); % Preallocate
for q = 1:Q
G(q,1) = 1i * 2*(a1(idx)/(omega(idx)*RI))* Z0d ...
* ( chi(q) / (chi(q)^2 - wk^2) )* ( besselj(l, wk*RI) / dbesselj(l, wk*RI) ); % this one is my calculated
end
% %G_qn
for q = 1:Q
for n = 2:N
G(q, n) = 1i * 2*(a1(idx)/(omega(idx)*RI))* Znd(n) ...
* ( chi(q) / (k(n)^2 + chi(q)^2) )* ( besseli(l, k(n)*RI) / dbesseli(l, k(n)*RI) ); % this one is my calculated
end
end
% % % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% %Define the right hand side vector
U = zeros(2*M + 2*N + Q, 1); % Full vector, M+N+M+N+Q elements
% % Block 1: U (size M = 4)
for i = 1:M
if i == 1
U(i, 1) = (sinh(wk*b))/((b*wk)*cosh(wk*h))*besselj(l, wk*RE) ; % Z_0^l
else
U(i, 1) = (2 *wk*(-1)^(i-1)*sinh(wk*b))/(b *(wk^2 +m(i)^2)*cosh(wk*h))*besselj(l, wk*RE); %Z_m^l
end
end
% % Block 2: Y (size N = 3)
for j = 1:N
if j == 1
U(j + M, 1) = -dbesselj(l, wk*RE) * (2*wk*h + sinh(2*wk*h)) /(cosh(wk*h)^2); % Y_0^l
else
U(j + M, 1) = 0; % Y_n^l
end
end
% % Block 3: X (size M = 4)
for i = 1:M
U(i + M + N, 1) = 0; % X_0^l, X_m^l
end
% Block 4: W (size N = 3)
for j = 1:N
U(j + M + N + M, 1) = 0; % W_0^l, W_n^l
end
for q = 1:Q
U(2*M + 2*N + q) = -(1i*b1/omega(idx)) * (2/(RI^2 * dbesselj(l, chi(q)*RI)));
end
% % Identity matrices
I_M = eye(M);
I_N = eye(N);
% Zero matrices
ZMM = zeros(M, M);
ZMN = zeros(M, N);
ZNM = zeros(N, M);
ZNN = zeros(N, N);
ZMQ = zeros(M,Q);
ZNQ = zeros(N,Q);
ZQM = zeros(Q,M);
ZQN = zeros(Q,N);
% Construct the full block matrix S
S = [ I_M, -A, ZMM, ZMN, ZMQ;
-B, I_N, -C, ZNN, ZNQ;
ZMM, ZMN, I_M, -D, ZMQ;
-E, ZNN, -F, I_N, -W;
ZQM, ZQN, ZQM, -G, H ];
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
solution_method = 0;
if solution_method == 0
if idx==1
psi0 = zeros(Q,1);
else
if exitflag == 1
psi0 = psi;
else
psi0 = zeros(Q,1);
end
end
f = @(psi)fun(psi,M,N,Q,S,U,Z0d,wk,RI,l,Znd,k,chi);
[psi,~,exitflag] = fsolve(f,psi0,optimset('Display','none'));
[~,d_vec,e_vec] = f(psi);
elseif solution_method == 1
converged = false;
psi = zeros(Q,1);
T_old = []; % <-- track the full solution from prev. iter
% T_old = zeros(2*M + 2*N + Q, 1); % previous unknowns (seed)
for iter = 1:max_iter
% (A) overwrite ONLY Block 5 of U from current psi
for q = 1:Q
U(2*M + 2*N + q) = -(1i*b1/omega(idx)) * (2/(RI^2 * dbesselj(l, chi(q)*RI))) * psi(q);
end
T = S \ U; % Solve the linear system
% T = pinv(S) * U; % Use pseudoinverse for stability
b_vec = T(1:M); % Coefficients b^l
a_vec = T(M+1 : M+N); % Coefficients a^l
c_vec = T(M+N+1 : 2*M+N); % Coefficients c^l
d_vec = T(2*M+N+1:2*M+2*N); % Coefficients d^l
e_vec = T(2*M+2*N+1:end); % (Q×1)
% (D) update psi for NEXT iteration from CURRENT coefficients
for q = 1:Q
integrand = @(r) abs(v_D(N,Q,r,Z0d,wk,RI,l,d_vec,Znd,k,e_vec,chi)) ...
.* v_D(N,Q,r,Z0d,wk,RI,l,d_vec,Znd,k,e_vec,chi) ...
.* besselj(l, chi(q)*r) .* r;
psi(q) = integral(integrand, 0, RI, 'AbsTol',1e-8, 'RelTol',1e-6);
end
% % % drawnow; % optional: flush output each iter
% === compact per-iteration print, only for T in [5,35] ===
% Tcur = 2*pi/omega(idx);
% if Tcur >= 5 && Tcur <= 35
% fprintf('idx %3d | iter %2d | ||psi|| = %.3e\n', idx, iter, norm(psi));
% drawnow;
% end
% % % % (E) convergence on ALL unknowns (T vs previous T)
% --- per-iteration diagnostics (optional) ---
if ~isempty(T_old)
diff_T = max(abs(T - T_old));
% Tcur = 2*pi/omega(idx);
% if Tcur >= 5 && Tcur <= 35
% fprintf('k0 idx %3d | iter %2d: max|ΔT| = %.3e, ||psi|| = %.3e\n', ...
% idx, iter, diff_T, norm(psi));
% drawnow; % ← add this to flush each iteration line too
% end
if diff_T < tol
converged = true;
break
end
end
T_old = T;
end % <<< end of iter loop
end
if solution_method==0
exitflag
elseif solution_method==1
converged
end
%T
% % ---- ONE summary line per frequency (place it here) ----
% iters_used(idx) = iter; % how many iterations this frequency needed
% fprintf('idx %3d (T=%6.3f s): iters=%2d\n', idx, 2*pi/omega(idx), iter);
% drawnow;
% -------------- end nonlinear iteration --------------
% % %%%%%%%%% Wave motion%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
term1 = (d_vec(1)*cosh(wk*h)^2) / ((2*wk*h + sinh(2*wk*h)) * dbesselj(l, wk*RI));
sum1 = 0;
for n =2:N
sum1 = sum1 + d_vec(n)*cos(k(n)*h)^2/ ((2*k(n)*h + sin(2*k(n)*h))* dbesseli(l, k(n)*RI));
end
phiUell = 0;
for q =1:Q
phiUell = phiUell +e_vec(q)* S_q(q);
end
eta0(idx) = abs(term1+sum1+phiUell);
end
% % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% % % Plotting Data
% % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% %
% %
% figure(2); clf;
% >>> CHANGES START (3): one figure = two taus + Dr. Chen only
%if old
%plot(2*pi./omega, eta0, 'k', 'LineWidth', 1.5); %%% this is T which is T= 2*pi/omega plotting against T
%else
plot(omega.^2 * RE / g, eta0, 'k', 'LineWidth', 1.5); %%% this is T which is T= 2*pi/omega plotting against T
%end
hold on; % Add extra plots on same figure
% % % %
% % % % % Now plot the 3 CSV experimental points
% scatter(data_chen(:,1), data_chen(:,2), 30, 'r', 'o', 'filled');
% % scatter(data_liu(:,1), data_liu(:,2), 30, 'g', 's', 'filled');
% scatter(data_exp(:,1), data_exp(:,2), 30, 'b', '^', 'filled');
xlabel('$T$', 'Interpreter', 'latex');
ylabel('$|\eta / (iA)|$', 'Interpreter', 'latex');
title('Wave motion amplitude for $R_E = 138$', 'Interpreter', 'latex');
legend({'$\tau=0.2$','Model test'}, ...
'Interpreter','latex','Location','northwest');
grid on;
% xlim([5 35]); % Match the reference plot
% ylim([0 4.5]); % Optional: based on expected peak height
%xlim([0 4]); % Match the reference plot
% ylim([0 7]); % Optional: based on expected peak height
% === plotting section ===
% diary off
elapsedTime = toc;
disp(['Time consuming = ', num2str(elapsedTime), ' s']);
function [res,d_vec,e_vec] = fun(psi,M,N,Q,S,U,Z0d,wk,RI,l,Znd,k,chi)
for q = 1:Q
U(2*M + 2*N + q) = U(2*M + 2*N + q) * psi(q);
end
T = S\U;
d_vec = T(2*M+N+1:2*M+2*N); % Coefficients d^l
e_vec = T(2*M+2*N+1:end); % (Q×1)
for q = 1:Q
integrand = @(r) abs(v_D(N,Q,r,Z0d,wk,RI,l,d_vec,Znd,k,e_vec,chi)) ...
.* v_D(N,Q,r,Z0d,wk,RI,l,d_vec,Znd,k,e_vec,chi) ...
.* besselj(l, chi(q)*r) .* r;
res(q) = psi(q) - integral(integrand, 0, RI, 'AbsTol',1e-8, 'RelTol',1e-6);
end
end
function out = dbesselk(l, z)
%DBESSELK Derivative of the modified Bessel function of the second kind
% out = dbesselk(l, z)
% Returns d/dz [K_l(z)] using the recurrence formula:
% K_l'(z) = -0.5 * (K_{l-1}(z) + K_{l+1}(z))
out = -0.5 * (besselk(l-1, z) + besselk(l+1, z));
end
function out = dbesselj(l, z)
%DBESSELJ Derivative of the Bessel function of the first kind
% out = dbesselj(l, z)
% Returns d/dz [J_l(z)] using the recurrence formula:
% J_l'(z) = 0.5 * (J_{l-1}(z) - J_{l+1}(z))
out = 0.5 * (besselj(l-1, z) - besselj(l+1, z));
end
function out = dbesseli(l, z)
%DBESSELI Derivative of the modified Bessel function of the first kind
% out = dbesseli(l, z)
% Returns d/dz [I_l(z)] using the recurrence formula:
% I_l'(z) = 0.5 * (I_{l-1}(z) + I_{l+1}(z))
out = 0.5 * (besseli(l-1, z) + besseli(l+1, z));
end
function out = dbesselh(l, z)
%DBESSELH Derivative of the Hankel function of the first kind
% out = dbesselh(l, z)
% Returns d/dz [H_l^{(1)}(z)] using the recurrence formula:
% H_l^{(1)'}(z) = 0.5 * (H_{l-1}^{(1)}(z) - H_{l+1}^{(1)}(z))
out = 0.5 * (besselh(l-1, 1, z) - besselh(l+1, 1, z));
end
function x = bessel0j(l,q,opt)
% a row vector of the first q roots of bessel function Jl(x), integer order.
% if opt = 'd', row vector of the first q roots of dJl(x)/dx, integer order.
% if opt is not provided, the default is zeros of Jl(x).
% all roots are positive, except when l=0,
% x=0 is included as a root of dJ0(x)/dx (standard convention).
%
% starting point for for zeros of Jl was borrowed from Cleve Moler,
% but the starting points for both Jl and Jl' can be found in
% Abramowitz and Stegun 9.5.12, 9.5.13.
%
% David Goodmanson
%
% x = bessel0j(l,q,opt)
k = 1:q;
if nargin==3 && opt=='d'
beta = (k + l/2 - 3/4)*pi;
mu = 4*l^2;
x = beta - (mu+3)./(8*beta) - 4*(7*mu^2+82*mu-9)./(3*(8*beta).^3);
for j=1:8
xnew = x - besseljd(l,x)./ ...
(besselj(l,x).*((l^2./x.^2)-1) -besseljd(l,x)./x);
x = xnew;
end
if l==0
x(1) = 0; % correct a small numerical difference from 0
end
else
beta = (k + l/2 - 1/4)*pi;
mu = 4*l^2;
x = beta - (mu-1)./(8*beta) - 4*(mu-1)*(7*mu-31)./(3*(8*beta).^3);
for j=1:8
xnew = x - besselj(l,x)./besseljd(l,x);
x = xnew;
end
end
end
% --- Local helper function for derivative of Bessel function ---
function dJ = besseljd(l, x)
dJ = 0.5 * (besselj(l - 1, x) - besselj(l + 1, x));
end
function v_D_val = v_D(N, Q, r, Z0d, wk, RI, l, d_vec, Znd, k, e_vec, chi)
% n = 0 mode
term1 = (d_vec(1) * Z0d * besselj(l, wk*r)) / (dbesselj(l, wk*RI));
% n >= 2 evanescent sum
sum1 = 0;
for nidx = 2:N
sum1 = sum1 + d_vec(nidx) * Znd(nidx) * besseli(l, k(nidx)*r) ...
/( dbesseli(l, k(nidx)*RI));
end
% q = 1..Q radial sum
sum2 = 0;
for qidx = 1:Q
sum2 = sum2 + e_vec(qidx) *chi(qidx)* besselj(l, chi(qidx)*r) / dbesselj(l, chi(qidx)*RI);
end
v_D_val = term1 + sum1 + sum2;
end
  38 Commenti
Javeria
Javeria il 25 Dic 2025 alle 13:20
@Torstenlicense('test','Distrib_Computing_Toolbox')
ans =
1
i tested my MATLAB
Torsten
Torsten il 25 Dic 2025 alle 14:02
Modificato: Torsten il 25 Dic 2025 alle 14:04
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.

Accedi per commentare.


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

Accedi per commentare.


idris
idris il 9 Gen 2026 alle 2:09
Modificato: idris il 9 Gen 2026 alle 2:16
Truly slow @Torsten, attached is the result
  6 Commenti
Javeria
Javeria il 12 Gen 2026 alle 2:39
@idris no, i did not solve my problem
Javeria
Javeria il 12 Gen 2026 alle 2:58
@Torsten see the below pdf and the equation number 24. They have the same non linear integral offcourse their expressions are changed but the non linear condition is same. They mentioned the method that how to solve it.

Accedi per commentare.


Javeria
Javeria il 12 Gen 2026 alle 4:17
@Torsten i have the following code for the numerical integration but i am not sure whether it will work for my case and then how i can use it for my case.
clear
% Consider f(x) = cos(x) * exp(sin(x)), i.e. d/dx [exp(sin(x))]
f = @(x) cos(x) .* exp(sin(x));
a = 0.6; b = 2.2;
% The analytical evaluation is exp(sin(b)) - exp(sin(a))
int_analytical = exp(sin(b)) - exp(sin(a));
% Set the number of grid points
N = 101;
% Choose a method 'T', 'S', or 'G'. Need N odd for 'S', but N can be odd or
% even for 'T' and 'G'
method = 'T';
[x,w] = get_quadrature_weights(a, b, N, method);
int_numerical = sum(w.*f(x));
disp(['Analytical evaluation: ' num2str(int_analytical, 10)])
Analytical evaluation: 0.4857117348
disp(['Numerical evaluation: ' num2str(int_numerical, 10)])
Numerical evaluation: 0.4856852319
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [x,w] = get_quadrature_weights(a, b, N, method)
% get quadrature nodes (x) and weights (w) with N points over the
% interval a < x < b. The method can be 'T' (trapezium rule), 'S'
% (Simpson's rule) or 'G' (Gauss-Legendre quadrature).
% Trapezium rule: the spacing between points is equal
% Simpson's rule: N needs to be odd
if strcmp(method, 'T')
dx = (b-a)/(N-1);
x = linspace(a,b,N);
w = dx/2 * [1, 2*ones(1,N-2), 1];
elseif strcmp(method, 'S')
if mod(N,2) == 0
error('N must be odd to use Simpson''s rule')
return
end
dx = (b-a)/(N-1);
x = linspace(a,b,N);
n = (N-1)/2;
A = [4*ones(1,n-1); 2*ones(1,n-1)];
w = dx/3 * [1, A(:).', 4, 1];
elseif strcmp(method, 'G')
[x,w] = lgwt(N,a,b);
else
error('Method needs to be "T", "S" or "G"')
x = NaN; w = NaN;
end
end
% ------------------------------------------------------------------
function [x,w]=lgwt(N,a,b)
% lgwt.m
%
% This script is for computing definite integrals using Legendre-Gauss
% Quadrature. Computes the Legendre-Gauss nodes and weights on an interval
% [a,b] with truncation order N
%
% Suppose you have a continuous function f(x) which is defined on [a,b]
% which you can evaluate at any x in [a,b]. Simply evaluate it at all of
% the values contained in the x vector to obtain a vector f. Then compute
% the definite integral using sum(f.*w);
%
% Written by Greg von Winckel - 02/25/2004
N=N-1;
N1=N+1; N2=N+2;
xu=linspace(-1,1,N1)';
% Initial guess
y=cos((2*(0:N)'+1)*pi/(2*N+2))+(0.27/N1)*sin(pi*xu*N/N2);
% Legendre-Gauss Vandermonde Matrix
L=zeros(N1,N2);
% Derivative of LGVM
Lp=zeros(N1,N2);
% Compute the zeros of the N+1 Legendre Polynomial
% using the recursion relation and the Newton-Raphson method
y0=2;
% Iterate until new points are uniformly within epsilon of old points
while max(abs(y-y0))>eps
L(:,1)=1;
Lp(:,1)=0;
L(:,2)=y;
Lp(:,2)=1;
for k=2:N1
L(:,k+1)=( (2*k-1)*y.*L(:,k)-(k-1)*L(:,k-1) )/k;
end
Lp=(N2)*( L(:,N1)-y.*L(:,N2) )./(1-y.^2);
y0=y;
y=y0-L(:,N2)./Lp;
end
% Linear map from[-1,1] to [a,b]
x=(a*(1-y)+b*(1+y))/2;
% Compute the weights
w=(b-a)./((1-y.^2).*Lp.^2)*(N2/N1)^2;
end
  18 Commenti
Javeria
Javeria il 15 Gen 2026 alle 14:21
@Torsten Thanks for your kind suggestion i tested your solvers by increasing N to test the convergence of series turncation also but when i increased it the computationall time is increased.
Torsten
Torsten il 15 Gen 2026 alle 19:33
Modificato: Torsten il 15 Gen 2026 alle 21:35
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

Accedi per commentare.

Categorie

Scopri di più su Linear Algebra in Help Center e File Exchange

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by