how to evaluate the integral of the expression involving bessel functions.

11 次查看（过去 30 天）
Javeria
Javeria 2025-7-18
评论： Torsten 2025-7-21
I have the following expression which i need to find.
then how i can find the integration? here this k is vector which contain k0,k1,k2, and same d and e is also a vector which contain d1,d2,d3 and e1, e2 ,e3 . and is also given. Then how i can evaluate this integration.
  2 个评论
Torsten
Torsten 2025-7-18
编辑：Torsten 2025-7-18
Why is psi a function of r ?
Which index p is meant when you write nu_D(r)*abs(nu_D(r))*besselj(chi_p^l * r) r as the integrand ?
Do you have functions to compute the derivatives of besselj and besseli ?
Why does your vector k has only 3 elements instead of 4 for k0,k1,k2 and k3 ?
Javeria
Javeria 2025-7-18
@Torsten yes i know the functions that how we can compute derivatives of bessel functions here we have k0,k1,k2 belongs to N=3.

请先登录，再进行评论。

请先登录，再回答此问题。

回答（2 个）

Shishir Reddy
Shishir Reddy 2025-7-18
编辑：Shishir Reddy 2025-7-18
Hi @Javeria
Kindly refer the following steps to numerically evaluate the given integral using MATLAB.
1. Define the Parameters:
% Given values
RI = 0.21;
k = [0.281, 0.2, 0.1];
d = [0.1, 0.3, 0.4]; % d0, d1, d2, d3 → d(1), d(2), d(3), d(4)
e = [0.3, 0.09, 0.7];
chi = [0.01, 0.03, 0.02];
2. Define ν_D(r):
nu_D = @(r) ...
d(1) * besselj(0, k(1)*r) / besselj(0, k(1)*RI) + ...
d(2) * besseli(1, k(2)*r) / besseli(1, k(2)*RI) + ...
d(3) * besseli(1, k(3)*r) / besseli(1, k(3)*RI) + ...
e(1) * besselj(1, chi(1)*r) / besselj(1, chi(1)*RI) + ...
e(2) * besselj(1, chi(2)*r) / besselj(1, chi(2)*RI) + ...
e(3) * besselj(1, chi(3)*r) / besselj(1, chi(3)*RI);
3. Define the Integrand Function:
chi1 = chi(1);
integrand = @(r) abs(nu_D(r)).^2 .* besselj(0, chi1 * r);
4. Use 'integral' to compute the result:
psi = integral(integrand, 0, RI);
disp(['ψ = ', num2str(psi)]);
For more information regarding 'integral' function, kindly refer the following documentation -
https://www.mathworks.com/help/matlab/ref/integral.html
I hope this helps.
  1 个评论
Javeria
Javeria 2025-7-18
@Shishir Reddy hello thanks for your response i might ask this question that how we can write this v_D(r) the way you write i understand but if we have a large N. then it's very difficult to write it.

请先登录，再进行评论。

Torsten
Torsten 2025-7-18
编辑：Torsten 2025-7-18
Ran in:
l = 0;
R_I = 2;
N = 3;
P = 3;
k = [0.28,0.2,0.1,0.05]; %[k0,k1,k2,k3]
e = [0.3,0.9,0.7]; %[e1,e2,e3]
d = [0.1,0.3,0.4,0.5]; %[d0,d1,d2,d3]
chi = [0.01,0.03,0.02]; %[chi1,chi2,chi3]
format long
psi = integral(@(r)fun(r,l,R_I,N,P,k,e,d,chi),0,R_I)
psi =
3.506800645644731e+02
function values = fun(r,l,R_I,N,P,k,e,d,chi)
nu_D_values = nu_D(r,l,R_I,N,P,k,e,d,chi);
chip = chi(end); %-> which p-index from p has to be taken here ?
values = nu_D_values.*abs(nu_D_values).*besselj(l,chip*r).*r;
end
function nu_D_values = nu_D(r,l,R_I,N,P,k,e,d,chi)
nu_D_values = d(1)*besselj(l,k(1)*r)/dbesselj(l,k(1)*R_I);
for n = 1:N
nu_D_values = nu_D_values + d(n+1)*besseli(l,k(n+1)*r)/dbesseli(l,k(n+1)*R_I);
end
for p = 1:P
nu_D_values = nu_D_values + e(p)*chi(p)*besselj(l,chi(p)*r)/dbesselj(l,chi(p)*R_I);
end
end
function value = dbesseli(nu,z)
value = 0.5*(besseli(nu-1,z) + besseli(nu+1,z));
%value = nu./z.*besseli(nu,z) + besseli(nu+1,z);
end
function value = dbesselj(nu,z)
value = 0.5*(besselj(nu-1,z) - besselj(nu+1,z));
%value = nu./z.*besselj(nu,z) - besselj(nu+1,z);
end
  5 个评论
显示 3更早的评论
Javeria
Javeria 2025-7-21
@Torsten Thanks for your response, I can share my full code. Actually i want to solve my system of equations for unknown coefficients b0,bm,a0,an,c0,cm,d0,dn,eq. However total i have 9 equations out of which 1 have non linear term due to this psi term. So i treat my this non linear term iteratively to solve my system of equations. However when i run my code i get zero vector for my unknowns which is physically not correct.
In my code i want to treat my psi iteratively that in 1st iteration it take the initial guess which i set to zero. And then onward it is updated by this expression. Below is my main script where i am calling helper functions for derivative of bessel functions and roots of bessel functions i set l=0 for the bessel functions mode.
The way i did for this psi is correct or not in my code?
% --- Main Script: Solve Non-Linear Moonpool Problem Iteratively ---
% Clear workspace, close all figures, and clear command window
close all;
clear all;
clc;
% --- 1. Define Global Constants and Parameters ---
% These define the truncation limits for the series summations in the equations.
N = 4; % Number of terms for 'n' modes (e.g., in Region 1/sum for 'a_n')
M = 4; % Number of terms for 'm' modes (e.g., in Region 2/sum for 'b_m')
Q = 3; % Number of terms for 'q' modes (e.g., in Region 3/sum for 'e_q')
% --- Physical Parameters ---
% These values are taken from the reference "Sphaier et al. (2007)".
RE = 47.5; % External radius of the main hull (m)
RI = 34.5; % Internal radius (m) - openness at the bottom of the moonpool
h = 200; % Water depth (m)
d = 38; % Draft (m)
g = 9.81; % Acceleration due to gravity (m/s^2)
b = (h-d); % 'b' as (h-d)
% --- Additional Constants ---
% These are constants specific to problem formulation.
gamma = 1.0; % empirical discharge coefficient term is in the range [0.5, 1.0]
tau = 0.2; % porosity ratio
% --- Bessel Function Order ---
% 'l' is the order of the Bessel functions used throughout your system of equations.
l = 0; % You can change this value as desired/required by the problem.
% This constant is derived from the formula (1 - tau) / (2 * gamma * tau^2)
b1 = (1 - tau) / (2 * gamma * tau^2);
% --- Iterative Solver Parameters ---
max_iter = 100;
tol = 1e-4;
converged = false;
% --- Wave Number and Group Velocity Calculation ---
X_kRE = linspace(5, 35, 100); % X-axis (kRE) matching the plot
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
I_given_wavenumber = 1; % Example flag, as in your code
% Loop through each dimensionless wavenumber (k0*RE) value
% For each X_kRE(i), we calculate the corresponding k0 (wk), omega, Cg, and
% the evanescent modes k(n)
for i = 1:length(X_kRE)
% Calculate the progressive wave mode wavenumber (k0),
%if I_given_wavenumber == 1
% wk = k0 = (k0*RE) / RE
wk = X_kRE(i) / RE; % wavenumber
omega = sqrt(g * wk * tanh(wk * h)); % wave frequency
Cg = (g*tanh(wk*h) + g*wk*h*(sech(wk*h))^2)*omega/(2*g*wk*tanh(wk*h)); %group velocity
a1 = 0.93 * (1 - tau) * Cg; %
fun_alpha = @(x) omega^2 + g * x * tan(x * h); % dispersion eqn
% end
for n = 1:N
if n == 1
k(n) = -1i * wk;
else
x0_n = [(2*n - 3)*pi/2/h, (2*n - 1)*pi/2/h];
k(n) = fzero(fun_alpha, mean(x0_n));
end
end
%%%%%%%%%%%%% Find derivative of Z0 at z=-d%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Z0d = (cosh(wk*h) * wk * sinh(wk*b))/(2*wk*h + sinh(2*wk*h));
%
%%%%%%%%%%%%% Find derivative of Zn at z=-d%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Znd = zeros(1, N); % or however many n you have
for n = 2:N % Typically n starts from 2 for these modes
Znd(n) = -k(n) * (cos(k(n)*h) * sin(k(n)*b)) / (2*k(n)*h + sin(2*k(n)*h));
end
for m=2:M
lambda(m) = ((m-1) * pi) / b;
end
% % Predefine identity and zero submatrices
O1 = zeros(M, M); % O1 == O6, size M×M
O2 = zeros(M, N); % O2 == O7, size M×N
O3 = zeros(M, Q); % O3 == O8 size M×Q
O4 = zeros(N, N); % O4 == O9, size N×N
O5 = zeros(N, Q); % size N×Q
O10 = zeros(Q, M); % O10 == O12, size Q×M
O11 = O5'; % O11 is the transpose of O5
I_M = eye(M); % size M×M
I_N = eye(N); % size N×N
%%%%%%%%%%%%%%%%%%%% Set Matrix A %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
A = zeros(M,N);
% For first element A00
A(1,1) = (cosh(wk*h) * sinh(wk*b)) / (wk*b * (2.*wk*h + sinh(2.*wk*h))) ...
* (besselh(l, 1, wk*RE) / wk*dbesselh(l, wk*RE));
% For 2nd Element A0n
for n = 2:N
A(1,n) = (cos(k(n)*h) * sin(k(n)*b)) / (k(n)*b * (2.*k(n)*h + sin(2.*k(n)*h)))...
* (besselk(l, k(n)*RE) /k(n)* dbesselk(l, k(n)*RE));
end
% for element Am0
for m = 2:M
A(m,1) = 2. * (cosh(wk*h) * (-1)^(m-1) * wk * sinh(wk*b)) / b*((wk^2 + (lambda(m))^2) ...
* (2.*wk*h + sinh(2.*wk*h))) * (besselh(l, 1, wk*RE) / wk*dbesselh(l, wk*RE));
end
%%%%%%%% for element Amn
for m = 2:M
for n = 2:N
A(m,n) = 2.* (cos(k(n)*h) * (-1)^(m-1)* k(n) * sin(k(n)*b)) / b*((k(n)^2 - (lambda(m))^2)...
* (2.*k(n)*h + sin(2.*k(n)*h))) * (besselk(l, k(n)*RE) / k(n)*dbesselk(l, k(n)*RE));
end
end
%%%%%%%%%%%%%%%%%% Matrix B %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
B = zeros(N, M); % Or whatever size you require
%%%%% Define B0%%%%%
B(1,1) = (4. * sinh(wk * b)) / (RE * log(RE / RI) * cosh(wk* h));
%%%%%%%% define Bom%%%%
for m = 2:M
Rml_prime_RE = lambda(m) * (besselk(l, lambda(m)*RI) * dbesseli(l, lambda(m)*RE) ...
- besseli(l, lambda(m)*RI) * dbesselk(l, lambda(m)*RE))...
/(besselk(l, lambda(m)*RI) * besseli(l, lambda(m)*RE) - besselk(l, lambda(m)*RE) * besseli(l, lambda(m)*RI));
B(1,m) = Rml_prime_RE * (4 * wk^2 * (-1)^(m-1) * sinh(wk*b)) / (cosh(wk*h) * (wk^2 + (lambda(m))^2));
end
%%%% define Bn0%%%%
for n = 2:N
B(n,1) = (4* sin(k(n)*b)) / (RE * log(RE/RI) * cos(k(n)*h));
end
%%%%%% Define Bnm %%%%%%%%
for n = 2:N
for m = 2:M
Rml_prime_RE = lambda(m) * (besselk(l, lambda(m)*RI) * dbesseli(l, lambda(m)*RE) ...
- besseli(l, lambda(m)*RI) * dbesselk(l, lambda(m)*RE))...
/(besselk(l, lambda(m)*RI) * besseli(l, lambda(m)*RE) - besselk(l, lambda(m)*RE) * besseli(l, lambda(m)*RI));
B(n,m) = Rml_prime_RE * (4* (k(n))^2 * (-1)^(m-1) * sin(k(n)*b)) / (cos(k(n)*h) * ((k(n))^2 - (lambda(m))^2));
end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Define Matrix C %%%%%%%%%%%%%%%%%%%%%%%%
C = zeros(N,M);
C(1,1) = -4 * sinh(wk*b) / (RE * log(RE/RI) * cosh(wk*h));
% %%%%% DEfine C0m%%%%%%
for m = 2:M
Rml_prime_star_RE = lambda(m) * (besseli(l, lambda(m)*RE) * dbesselk(l, lambda(m)*RE) ...
- besselk(l, lambda(m)*RE) * dbesseli(l, lambda(m)*RE))...
/(besselk(l, lambda(m)*RI) * besseli(l, lambda(m)*RE) - besselk(l, lambda(m)*RE) * besseli(l, lambda(m)*RI));
C(1,m) = Rml_prime_star_RE * (4 * wk^2 * (-1)^(m-1) * sinh(wk*b)) / (cosh(wk*h) * (wk^2 + (lambda(m))^2));
end
%%%%%%% Define Cn0%%%%%%
for n = 2:N
C(n,1) = -4 * sin(k(n)*b) / (RE * log(RE/RI) * cos(k(n)*h));
end
%%%%%% Define Cnm%%%%%%
for n = 2:N
for m =2:M
Rml_prime_star_RE = lambda(m) * (besseli(l, lambda(m)*RE) * dbesselk(l, lambda(m)*RE) ...
- besselk(l, lambda(m)*RE) * dbesseli(l, lambda(m)*RE))...
/(besselk(l, lambda(m)*RI) * besseli(l, lambda(m)*RE) - besselk(l, lambda(m)*RE) * besseli(l, lambda(m)*RI));
C(n,m) = Rml_prime_star_RE * (4 * (k(n))^2 * (-1)^(m-1) * sin(k(n)*b)) / (cos(k(n)*h) * ((k(n))^2 - (lambda(m))^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) / wk*dbesselj(l, wk*RI));
%%%%% write D0n %%%%%%%%%%
for n =2:N
D(1,n) = (cos(k(n)*h) * sin(k(n)*b)) / (k(n)*b * (2.*k(n)*h + sin(2.*k(n)*h))) * (besseli(l, k(n)*RI) / k(n)*dbesseli(l, k(n)*RI));
end
%%%%%% write Dm0 %%%%%%%%%
for m = 2:M
D(m,1) = (2 * cosh(wk*h) * (-1)^(m-1) * wk * sinh(wk*b)) /(b * (2.*wk*h + sinh(2.*wk*h)) * (wk^2 + (lambda(m))^2))...
*(besselj(l, wk*RI) /wk* dbesselj(l, wk*RI));
end
%%%%% Define Dmn%%%%%%%%
for m = 2:M
for n = 2:N
D(m,n) = (2 * cos(k(n)*h) * (-1)^(m-1) * k(n) * sin(k(n)*b)) /(b * (2.*k(n)*h + sin(2.*k(n)*h)) * ((k(n))^2 - (lambda(m))^2))...
*(besseli(l, k(n)*RI) /k(n)* dbesseli(l, k(n)*RI));
end
end
%%%%%%%%%%%%%%%%%%%%%%%%%% Define Matrix E %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
E = zeros(N, M); % Preallocate the matrix E
%%%%% Define E00%%%%%%%%%%%%
E(1,1) = (4 * sinh(wk*b)) / (RI * log(RE/RI) * cosh(wk*h));
%%%% Define Eom %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for m = 2:M
Rml_prime_RI = lambda(m)*(besselk(l, lambda(m)*RI) * dbesseli(l, lambda(m)*RI) ...
- besseli(l, lambda(m)*RI) * dbesselk(l, lambda(m)*RI))...
/ (besselk(l, lambda(m)*RI) * besseli(l, lambda(m)*RE) - besselk(l, lambda(m)*RE) * besseli(l, lambda(m)*RI));
E(1,m) = Rml_prime_RI * (4 * wk^2 * (-1)^(m-1) * sinh(wk*b)) / (cosh(wk*h) * (wk^2 + lambda(m)^2));
end
%%%%%% Define Eno%%%%%%%%%%%%%%
for n = 2:N
E(n,1) = (4 * sin(k(n)*b)) / (RI * log(RE/RI) * cos(k(n)*h));
end
for n = 2:N
for m = 2:M
Rml_prime_RI = lambda(m)*(besselk(l, lambda(m)*RI) * dbesseli(l, lambda(m)*RI) ...
- besseli(l, lambda(m)*RI) * dbesselk(l, lambda(m)*RI))...
/ (besselk(l, lambda(m)*RI) * besseli(l, lambda(m)*RE) - besselk(l, lambda(m)*RE) * besseli(l, lambda(m)*RI));
E(n,m) = Rml_prime_RI * (4 * k(n)^2 * (-1)^(m-1) * sin(k(n)*b)) / (cos(k(n)*h) * (k(n)^2 - lambda(m)^2));
end
end
%%%%%%%%%%%%%%%%%%%%%%%% Now write matrix F %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
F = zeros(N,M);
%%%%%%%% Define F00 %%%%%%%%%%
F(1,1) = (-4 * sinh(wk*b)) / (RI * log(RE/RI) * cosh(wk*h));
%%%%%% Define F0m%%%%
for m = 2:M
Rml_star_prime_RE = lambda(m)*(besseli(l, lambda(m)*RE) * dbesselk(l, lambda(m)*RI) ...
- besselk(l, lambda(m)*RE) * dbesseli(l, lambda(m)*RI))/((besselk(l, lambda(m)*RI) ...
* besseli(l, lambda(m)*RE) - besselk(l, lambda(m)*RE) * besseli(l, lambda(m)*RI)));
F(1,m) = Rml_star_prime_RE * (4 * wk^2 * (-1)^(m-1) * sinh(wk*b)) / (cosh(wk*h) * (wk^2 + lambda(m)^2));
end
%%%%% Defin Fn0%%%%%%
for n = 2:N
F(n,1) = (-4 * sin(k(n)*b)) / (RI * log(RE/RI) * cos(k(n)*h));
end
%%%%%%% Define Fnm %%%%%%%
for n = 2:N
for m = 2:M
Rml_star_prime_RE = lambda(m)*(besseli(l, lambda(m)*RE) * dbesselk(l, lambda(m)*RI) ...
- besselk(l, lambda(m)*RE) * dbesseli(l, lambda(m)*RI))/((besselk(l, lambda(m)*RI) ...
* besseli(l, lambda(m)*RE) - besselk(l, lambda(m)*RE) * besseli(l, lambda(m)*RI)));
F(n,m) = Rml_star_prime_RE * (4 * k(n)^2 * (-1)^(m-1) * sin(k(n)*b)) / (cos(k(n)*h) * (k(n)^2 - lambda(m)^2));
end
end
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Calling the bessel0j(l,q,option) to find the positive roots of Bessel functions
roots = bessel0j(l, Q); % [x_1^l, x_2^l, ..., x_Q^l] (1 x Q vector)
chi = roots / RI; % (1 x Q vector) -- only if you need chi_q^l
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%% Write Matrix W %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
W = zeros(N,Q);
%%W_{0q}
for q = 1:Q
W(1, q) = -4*wk*(chi(q)*cosh(wk*b)-wk*coth(chi(q)*b)*sinh(wk*b)) ...
/ ((chi(q)^2-wk^2) * cosh(wk*h));
end
%%% Write W_{nq}
for n = 2:N
for q = 1:Q
W(n, q) = -4*k(n)*( chi(q)*cos(k(n)*b) + k(n)*coth(chi(q)*b)*sin(k(n)*b)) ...
/ ( cos(k(n)*h)*(chi(q)^2 + k(n)^2) );
end
end
% %%%%%%%%%%%%%%%%%%%%%%%Find E1_{q}%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for q = 1:Q
f2 = omega^2 / g;
E1(q) = ( chi(q)*cosh(chi(q)*h) - f2*sinh(chi(q)*h) ) ...
/ ( f2*cosh(chi(q)*d) - chi(q)*sinh(chi(q)*d) );
end
G = zeros(Q, N); % Preallocate
for q = 1:Q
% Anonymous function for integrand
f1 = @(r) (besselj(l, wk*r) .* besselj(l, chi(q)*r) .* r);
% Numerical integration over [0, RI]
int1(q, 1) = integral(f1, 0, RI);
% G_{q0}^l formula (assuming Z0d is constant, dbesselj is your derivative of Jl)
G(q, 1) = (-1i * a1 / omega)*Z0d* int1(q, 1) / dbesselj(l, wk*RI);
end
%% G_{qn}
for q = 1:Q
for n = 2:N
% Anonymous function for integrand
f2 = @(r) (besseli(l, k(n)*r) .* besselj(l, chi(q)*r) .* r);
% Numerical integration over [0, RI]
int2(q, n) = integral(f2, 0, RI);
% G_{qn}^l formula (assuming Znd is constant, dbesselj is your derivative of Jl)
G(q, n) = (-1i * a1 / omega)*Znd(n) * int2(q, n) / dbesseli(l, k(n)*RI);
end
end
% find H_{q}
for q = 1:Q
% Compute the diagonal value for H at index q
H(q) = (2 / (sinh(2 * chi(q) * (h - d))) - E1(q) + 1i * a1 * chi(q) / omega ) ...
* RI^2*( besselj(l+1, chi(q) * RI)^2 )/(2 * dbesselj(l, chi(q) * RI));
end
H_I = diag(H); % Diagonal matrix
T_old = zeros(2*M + 2*N + Q, 1); % Preallocate for previous unknowns
% --- Initial guess: |v_D(r)| = 0, so psi = 0
psi = zeros(Q, 1);
for iter = 1:max_iter
%Define the right hand side vector
U = zeros(2*M + 2*N + Q, 1); % Full vector, M+N+M+N+Q elements
for m = 1:M
if m == 1
% Z_0^l term (first entry in Block 1)
U(m,1) = (besselj(0, wk*RE) * sinh(wk*b)) / (b*wk * cosh(wk*h));
else
% Z_m^l terms (m = 2,3,...,M)
U(m,1) = (2 * besselj(l, wk*RE) * sinh(wk*b) * wk * (-1)^(m-1)) ...
/ (b* cosh(wk*h)*(wk^2 + lambda(m)^2));
end
end
% Block 2: Y (size N = 3)
for n = 1:N
if n == 1
U(n + M, 1) = -wk*dbesselj(l, wk*RE) * (2*wk*h + sinh(2*wk*h)) /(cosh(wk*h)^2); % Y_0^l
else
U(n + M, 1) = 0; % Y_n^l
end
end
% Block 3: X (size M)
for m = 1:M
U(m + M + N, 1) = 0; % X_0^l, X_m^l
end
% Block 4: W (size N)
for n = 1:N
U(n + M + N + M, 1) = 0; % K_0^l, K_n^l
end
% Block 5: Nonlinear part
for q = 1:Q
U(2*M + 2*N + q) = -1i * b1 / omega * psi(q);
end
%
%%%%%%%%%% DEFINE THE MATRIX S%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
S = [ I_M, -A, O1, O2, O3;
-B, I_N, -C, O4, O5;
O1, O2, I_M, -D, O3;
-E, O4, -F, I_N, -W;
O10, O11, O10, -G, H_I];
% Solve for the unknowns vector T (which holds [bm; an; cm; dn; eq])
T = pinv(S) * U; % If S may be singular or not square
% === (C) Extract unknowns from T for current iteration ===
bm = T(1:M);
an = T(M+1:M+N);
cm = T(M+N+1:2*M+N);
dn = T(2*M+N+1:2*M+2*N);
eq = T(2*M+2*N+1:end);
% === (D) Update psi(q) for next iteration ===
for q = 1:Q
integrand = @(r) abs(v_D(N, Q, r, Z0d, wk, RI, l, dn, Znd, k, eq, chi)) ...
.* v_D(N, Q, r, Z0d, wk, RI, l, dn, Znd, k, eq, chi) ...
.* besselj(l, chi(q)*r) .* r;
psi(q) = integral(integrand, 0, RI, 'AbsTol', 1e-8, 'RelTol', 1e-8);
end
if iter > 1
diff_T = max(abs(T - T_old));
fprintf('Iteration %d: max(|T-T_old|) = %.3e\n', iter, diff_T);
if diff_T < tol
converged = true;
fprintf(' Solution converged after %d inner iterations.\n', iter);
break;
end
end
T_old = T;
end
% --- Store results for the current 'i' (X_kRE) case ---
% Use 'i' as the index for storing results in cell arrays.
if converged
bm_all{i} = bm;
an_all{i} = an;
cm_all{i} = cm;
dn_all{i} = dn;
eq_all{i} = eq;
T_all{i} = T;
psi_all{i}= psi;
else
warning('Did not converge for X_kRE point %d. Max iterations reached.', i);
bm_all{i} = NaN; an_all{i} = NaN; cm_all{i} = NaN; dn_all{i} = NaN; eq_all{i} = NaN;
T_all{i} = NaN;
psi_all{i}= NaN;
end
end
Below is the helper functions which i used for my main code.
For the positive roots of bessel functions
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
For the derivative of bessel function
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 = 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 = 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 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
Expression for v_D(r)
function v_D_val = v_D(N, Q, r, Z0d, wk, RI, l, dn, Znd, k, eq, chi)
% Computes v_D(r) for arbitrary r, given coefficients and parameters.
% Inputs:
% N - number of n modes
% Q - number of q modes
% r - radial coordinate (scalar)
% Z0, wk, RI, l, dn, Zn, k, eq, chi - see above documentation
% First term: n=0 mode
term1 = (dn(1) * Z0d * besselj(l, wk*r)) / (wk*dbesselj(l, wk*RI));
% Second term: sum over n=2:N
sum1=0;
for nidx = 2:N % For N=3, this runs for nidx = 2 and nidx = 3
sum1= sum1 + dn(nidx) * Znd(nidx) * besseli(l, k(nidx)*r) /(k(nidx)* dbesseli(l, k(nidx)*RI));
end
% Third term: sum over q=1:Q
sum2 =0;
for qidx = 1:Q
sum2= sum2+ eq(qidx) * chi(qidx) * besselj(l, chi(qidx)*r) / (chi(qidx)*dbesselj(l, chi(qidx)*RI));
end
% Sum all terms to get the final v_D(r) value
v_D_val = term1 + sum1 + sum2;
end
Torsten
Torsten 2025-7-21
The first thing that I'd do is to output the matrices A, B, C, D E, F, W, G and H_I and inspect where the entries in the order of 1e13 - 1e29 stem from in your computations.

请先登录，再进行评论。

类别

在 Help Center 和 File Exchange 中查找有关 Special Functions 的更多信息

标签

Community Treasure Hunt

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

Start Hunting!

Translated by