Integrations of Bessel functions
Show older comments
I have these two integrals and i want to find its analytical integral i tried in MATLAB symbolic tool box but i could not.
How to find this integral analytically.
Accepted Answer
For the first integral I_1, see
For this example, the closed-form expression works:
f = @(x)besselj(0,2*x).*besselj(0,3*x).*x;
integral(f,0,4,'AbsTol',1e-12,'RelTol',1e-12)
4/(2^2-3^2)*(3*besselj(0,2*4)*besselj(-1,3*4)-2*besselj(-1,2*4)*besselj(0,3*4))
But I can't reproduce the result claimed under "One obtains finally":
f = @(x)besseli(0,2*x).*besselj(0,3*x).*x;
integral(f,0,4,'AbsTol',1e-12,'RelTol',1e-12)
4*3/(2^2+3^2)*besselj(1,3*4)*besseli(0,2*4)
Wolframalpha gives a solution for I_2 for l = 0 and l = 1:
More Answers (1)
Umar
on 16 Aug 2025
0 votes
Hi @Javeria,
I read your comments and implemented the analytical formulas for the two Bessel function integrals directly in MATLAB. The script computes their exact values without using the symbolic toolbox, so it provides the analytical results you were looking for. This should fully answer your questions.
Implementation
%========================================================== % Bessel product integrals (analytic evaluation, no toolbox) %==========================================================
% Integral 1: % I = ∫_0^a r J0(nu r) I0(mu r) dr % General closed form: % I = (a/(nu^2+mu^2)) * [ nu*J1(nu a) I0(mu a) + mu*J0(nu a) I1(mu a) ] % % If J0(nu a) = 0 (typical boundary condition), it reduces to: % I = (a*nu/(nu^2+mu^2)) * J1(nu a) I0(mu a)
function val = besselJI0_int(nu, mu, a) term1 = nu * besselj(1, nu*a) * besseli(0, mu*a); term2 = mu * besselj(0, nu*a) * besseli(1, mu*a); val = a/(nu^2 + mu^2) * (term1 + term2); end
% ------------------------------------------------------------
% Integral 2:
% I = ∫_0^R r J_l(k0 r) I_l(chi r) dr
% General closed form:
% I = (R/(k0^2+chi^2)) * [ k0 J_{l+1}(k0 R) I_l(chi R)
% + chi J_l(k0 R) I_{l+1}(chi R) ]
function val = besselJI_int(l, k0, chi, R) term1 = k0 * besselj(l+1, k0*R) * besseli(l, chi*R); term2 = chi * besselj(l, k0*R) * besseli(l+1, chi*R); val = R/(k0^2 + chi^2) * (term1 + term2); end
% ============================================================ % Example usage % ============================================================
% Parameters nu = 2.5; % example value for nu mu = 1.2; % example value for mu a = 1.0; % upper limit
k0 = 3.0; % example k0 chi = 1.5; % example chi R = 1.0; % upper limit l = 0; % order
% Compute integrals I1 = besselJI0_int(nu, mu, a); I2 = besselJI_int(l, k0, chi, R);
% Display results
fprintf('Integral 1 result = %.6f\n', I1);
fprintf('Integral 2 result = %.6f\n', I2);
Results
Furthermore, My solution implements the exact analytical formulas for integrals of products of Bessel J and modified Bessel I functions, allowing direct computation without the symbolic toolbox. I also appreciate @Torsten’s comments—they correctly point out that formulas for J⋅J integrals don’t directly carry over to J⋅I integrals. His observations helped clarify why a naive application of the J⋅J formula wouldn’t match the numerical result. My script builds on this understanding to provide the correct analytical values for the J⋅I case.
17 Comments
@Umar @Torsten Thanks for the answer actually I have just a fortan code as a refrence for my problem and i just implement their analytical results but i did not get the desired results i have a dissipative boundary condition and i want to simplify it to incoporate it inside in my code for the unknown coefficeint c_n. We make use of the orthogonal bessel functions which will help us to take our the unknown equation.
where
and
in fortan they write this one as
do iM=1, nbM
ijM=mij(iM+2*(nbN+1))
Lm0=ZpHmD(0)*knR(0)*kmR(iM)*BJsBJp(0) !This one is for J_l(k0r).J_l(\gammar)
Lm0=Lm0/(knR(0)*knR(0)-kmR(iM)*kmR(iM))
Lm0=Lm0*2.0D0*dissip
ijM=ijM+1+nbN+1
matA(ijM)=ZI*Lm0
do iN=1, nbN
Lm0=ZpHmD(iN)*knR(iN)*kmR(iM)*BIsBIp(iN)
Lm0=Lm0/(knR(iN)*knR(iN)+kmR(iM)*kmR(iM))
Lm0=Lm0*2.0D0*dissip
ijM=ijM+1
matA(ijM)=ZI*Lm0
end do
Lm0=4.0D0*exp2kmHmD(iM)/(1.0D0-exp2kmHmD(iM)*exp2kmHmD(iM))-CCm(iM)
ijM=ijM+iM
matA(ijM)=Lm0+ZI*dissip*kmR(iM)
end do
The confusion i have with R^2/2 when for simplification we multiply both side by this one. When i do simplifications then i have left with 2/R term on right hand side 1st and 2nd term but in fortan they do not write this R. So i am a little bit confused about this.
Umar
on 17 Aug 2025
Hi @ Javeria,
I’m not very familiar with Fortran, but this link might help you see how similar code works or how to translate it to MATLAB:
Let me know if I can help with anything else.
Umar
on 17 Aug 2025
I might have pasted the wrong link, here is the correct one
Torsten
on 17 Aug 2025
Take a look at relation (2) under
Do you use it in the form as stated (with k_1 = mu_n/l and k_2 = mu_k/l) or with k_1 = mu_n and k_2 = mu_k ?
Can you show us your computations where you noticed the discrepancy about the R^2/2 ?
Javeria
on 17 Aug 2025
Umar
on 17 Aug 2025
Hi @Javeria,
Thanks for sharing your work. Looking at your derivation and the Brown reference Torsten mentioned, I can see where the R-squared over 2 discrepancy is coming from.
Regarding Torsten's question about the scaling - you should definitely use the scaled version where you divide the eigenvalue by the length, as stated in relation (2). In your handwritten work (Dissipative -integral.pdf), I notice you're using the bare eigenvalues in your Bessel functions, but you need to include the "r divided by R" scaling factor inside the Bessel function arguments.
For example, instead of writing:"Bessel function of the first kind of order zero, evaluated at chi-m times r", you should write: "Bessel function of the first kind of order zero, evaluated at chi-m times r divided by R"
Here's how to rework your derivation step by step:
Start fresh with proper scaling: Go through every single Bessel function in your equations and make sure they all have the "times r divided by R" form inside them.
Set up the series expansion correctly: Instead of jumping straight to orthogonality, first write your unknown function as an infinite sum of properly scaled Bessel functions with unknown coefficients.
Apply orthogonality systematically: Multiply both sides of your main equation by a properly scaled Bessel function, then integrate from zero to R. The orthogonality property will make most terms disappear, leaving you with just one coefficient to solve for.
Let the normalization appear naturally: Don't try to force the R-squared over 2 factor into your algebra. It will appear automatically when you evaluate the remaining integral on the left side.
Check dimensions throughout: Every term in your equation should have the same physical units. If you're getting extra R factors, it usually means the scaling isn't consistent somewhere.
Does this approach make sense for redoing the calculation?
I get the same results for the integrals as you do, except for the integral
integral_{r=0}^{r=R} r * J_l(gamma_m^l * r) * I_l(k_n*r) dr
which should be
-gamma_m^l * R/(k_n^2 + gamma_m^l^2) * J_l'(gamma_m^l*R) * I_l(k_n*R)
I guess (thus your result with a minus sign).
Concerning the remaining R: What are the physical units of Z_0' and Z_n' ? Is it 1/m ?
A tip from my side:
You can easily test the validity of analytical formulae for certain integrals by a numerical check (as done in the enclosed code for the above integral relation).
format long
m = 20; % maximum number of roots required
kind = 1; % We are dealing with bessel functions of the 1st kind
l0 = 0; % order of bessel function = 0
x = besselzero(l0,m,kind); % compute the m positive roots j_{0,1:20}
%besselj(0,x) % Check if J0(x) = 0
R = 8;
x = x/R;
kn = 4;
f = @(r)r.*besselj(l0,x(12)*r).*besseli(l0,kn*r);
integral(f,0,R,'RelTol',1e-12,'AbsTol',1e-12)
-R*x(12)/(kn^2+x(12)^2)*dbesselj(l0,x(12)*R)*besseli(l0,kn*R)
function x = besselzero(n, k, kind)
% besselzero calculates the zeros of Bessel function of the first and second kind
%
% x = besselzero(n)
% x = besselzero(n, k)
% x = besselzero(n, k, kind)
%
%% Inputs
% * *n* - The order of the bessel function. n can be a scalar, vector, or
% matrix. n can be positive, negative, fractional, or any
% combinaiton. abs(n) must be less than or equal to
% 146222.16674537213 or 370030.762407380 for first and second kind
% respectively. Above these values, this algorithm will not find the
% correct zeros because of the starting values therefore an error is
% thrown instead.
% * k - The number of postive zeros to calculate. When k is not supplied,
% k = 5 is the default. k must be a scalar.
% * kind - kind is either 1 or 2. When kind is not supplied, default is
% kind = 1.
%
%% Outputs
% * x - The calculated zeros. size(x) = [size(n) k].
%
%% Description
% besselzero calculates the first k positive zeros of nth order bessel
% function of the first or second kind. Note, that zero is not included as
% the first zero.
%
%% Algorithm
% the first three roots of any order bessel can be approximated by a simple
% equations. These equations were generated using a least squares fit of
% the roots from orders of n=0:10000. The approximation is used to start
% the iteration of Halley's method. The 4th and higher roots can be
% approximated by understanding the roots are regularly spaced for a given
% order. Once the 2nd and 3rd roots are found, the spacing can be
% approximated by the distance between the 2nd and 3rd root. Then again
% Halley's method can be applied to precisely locate the root.
%%
% Because the algorithm depends on good guesses of the first three zeros,
% if the guess is to far away then Halley's method will converge to the
% wrong zero which will subsequently cause any other zero to be incorrectly
% located. Therefore, a limit is put on abs(n) of 146222.16674537213 and
% 370030.762407380 for first and second kind respectively. If n is
% specified above these limits, then an error is thrown.
%
%% Example
% n = (1:2)';
% k = 10;
% kind = 1;
% z = besselzero(n, k, kind);
% x = linspace(0, z(end), 1000);
% y = nan(2, length(x));
% y(1,:) = besselj(n(1), x);
% y(2,:) = besselj(n(2), x);
% nz = nan(size(z));
% nz(1,:) = besselj(n(1), z(1,:));
% nz(2,:) = besselj(n(2), z(2,:));
% plot(x, y, z, nz,'kx')
%
% Originally written by
% Written by: Greg von Winckel - 01/25/05
% Contact: gregvw(at)chtm(dot)unm(dot)edu
%
% Modified, Improved, and Documented by
% Jason Nicholson 2014-Nov-06
% Contact: jashale@yahoo.com
%% Change Log
% * Original release. 2005-Jan-25, Greg von Winckel.
% * Updated Documentation and commented algorithm. Fixed bug in finding the
% the first zero of the bessel function of the second kind. Improved speed by
% factor of 20. 2014-Nov-06, Jason Nicholson.
%
%% Input checking
assert(nargin >=1 | nargin <=3,'Wrong number of input arguments.');
% Take care of default cases of k and kind
if nargin < 2
k = 5;
end
if nargin < 3
kind = 1;
end
assert(isscalar(kind) & any(kind == [1 2]), '''kind''must be a scalar with value 1 or 2 only.');
assert(isscalar(k) & fix(k)==k & k>0, 'k must a positive scalar integer.');
assert(all(isreal(n(:))), 'n must be a real number.');
% negative orders have the same roots as the positive orders
n = abs(n);
% Check for that n is less than the ORDER_MAX
if kind==1
ORDER_MAX = 146222.16674537213;
assert(all(n <= ORDER_MAX), 'all n values must be less than or equal %6.10f for kind=1.', ORDER_MAX);
elseif kind==2
ORDER_MAX = 370030.762407380;
assert(all(n(:) <= ORDER_MAX), 'all n values must be less than or equal %6.10f for kind=2.', ORDER_MAX);
end
%% Setup Arrays
% output size
nSize = size(n);
if nSize(end) ==1
outputSize = [nSize(1:end-1) k];
else
outputSize = [nSize k];
end
% number of orders for each kth root
nOrdersPerRoot = prod(outputSize(1:end-1));
x = nan(outputSize);
%% Solve for Roots
switch kind
case 1
% coefficients and exponent are from least squares fitting the k=1,
% n=0:10000.
coefficients1j = [0.411557013144507;0.999986723293410;0.698028985524484;1.06977507291468];
exponent1j = [0.335300369843979,0.339671493811664];
% guess for k = 1
x((1:nOrdersPerRoot)') = coefficients1j(1) + coefficients1j(2)*n(:) + coefficients1j(3)*(n(:)+1).^(exponent1j(1)) + coefficients1j(4)*(n(:)+1).^(exponent1j(2));
% find first zero
x((1:nOrdersPerRoot)') = arrayfun(@(n, x0) findzero(n, 1, x0, kind), n(:), x((1:nOrdersPerRoot)'));
if k >= 2
% coefficients and exponent are from least squares fitting the k=2,
% n=0:10000.
coefficients2j = [1.93395115137444;1.00007656297072;-0.805720018377132;3.38764629174694];
exponent2j = [0.456215294517928,0.388380341189200];
% guess for k = 2
x((nOrdersPerRoot+1:2*nOrdersPerRoot)') = coefficients2j(1) + coefficients2j(2)*n(:) + coefficients2j(3)*(n(:)+1).^(exponent2j(1)) + coefficients2j(4)*(n(:)+1).^(exponent2j(2));
% find second zero
x((nOrdersPerRoot+1:2*nOrdersPerRoot)') = arrayfun(@(n, x0) findzero(n, 2, x0, kind), n(:), x((nOrdersPerRoot+1:2*nOrdersPerRoot)'));
end
if k >= 3
% coefficients and exponent are from least squares fitting the k=3,
% n=0:10000.
coefficients3j = [5.40770803992613;1.00093850589418;2.66926179799040;-0.174925559314932];
exponent3j = [0.429702214054531,0.633480051735955];
% guess for k = 3
x((2*nOrdersPerRoot+1:3*nOrdersPerRoot)') = coefficients3j(1) + coefficients3j(2)*n(:) + coefficients3j(3)*(n(:)+1).^(exponent3j(1)) + coefficients3j(4)*(n(:)+1).^(exponent3j(2));
% find second zero
x((2*nOrdersPerRoot+1:3*nOrdersPerRoot)') = arrayfun(@(n, x0) findzero(n, 3, x0, kind), n(:), x((2*nOrdersPerRoot+1:3*nOrdersPerRoot)'));
end
case 2
% coefficients and exponent are from least squares fitting the k=1,
% n=0:10000.
coefficients1y = [0.0795046982450635;0.999998378297752;0.890380645613825;0.0270604048106402];
exponent1y = [0.335377217953294,0.308720059086699];
% guess for k = 1
x((1:nOrdersPerRoot)') = coefficients1y(1) + coefficients1y(2)*n(:) + coefficients1y(3)*(n(:)+1).^(exponent1y(1)) + coefficients1y(4)*(n(:)+1).^(exponent1y(2));
% find first zero
x((1:nOrdersPerRoot)') = arrayfun(@(n, x0) findzero(n, 1, x0, kind), n(:), x((1:nOrdersPerRoot)'));
if k >= 2
% coefficients and exponent are from least squares fitting the k=2,
% n=0:10000.
coefficients2y = [1.04502538172394;1.00002054874161;-0.437921325402985;2.70113114990400];
exponent2y = [0.434823025111322,0.366245194174671];
% guess for k = 2
x((nOrdersPerRoot+1:2*nOrdersPerRoot)') = coefficients2y(1) + coefficients2y(2)*n(:) + coefficients2y(3)*(n(:)+1).^(exponent2y(1)) + coefficients2y(4)*(n(:)+1).^(exponent2y(2));
% find second zero
x((nOrdersPerRoot+1:2*nOrdersPerRoot)') = arrayfun(@(n, x0) findzero(n, 2, x0, kind), n(:), x((nOrdersPerRoot+1:2*nOrdersPerRoot)'));
end
if k >= 3
% coefficients and exponent are from least squares fitting the k=3,
% n=0:10000.
coefficients3y = [3.72777931751914;1.00035294977757;2.68566718444899;-0.112980454967090];
exponent3y = [0.398247585896959,0.604770035236606];
% guess for k = 3
x((2*nOrdersPerRoot+1:3*nOrdersPerRoot)') = coefficients3y(1) + coefficients3y(2)*n(:) + coefficients3y(3)*(n(:)+1).^(exponent3y(1)) + coefficients3y(4)*(n(:)+1).^(exponent3y(2));
% find second zero
x((2*nOrdersPerRoot+1:3*nOrdersPerRoot)') = arrayfun(@(n, x0) findzero(n, 3, x0, kind), n(:), x((2*nOrdersPerRoot+1:3*nOrdersPerRoot)'));
end
otherwise
error('Code should never get here.');
end
if k >= 4
for iRoot = 4:k
% guesses for remaining roots x[k] = rootSpacing + x[k-1]
x(((iRoot-1)*nOrdersPerRoot+1:iRoot*nOrdersPerRoot)') = 2*x(((iRoot-2)*nOrdersPerRoot+1:(iRoot-1)*nOrdersPerRoot)')- x(((iRoot-3)*nOrdersPerRoot+1:(iRoot-2)*nOrdersPerRoot)');
% find the remaining zeros
x(((iRoot-1)*nOrdersPerRoot+1:iRoot*nOrdersPerRoot)') = arrayfun(@(n, x0) findzero(n, k, x0, kind), n, x(((iRoot-1)*nOrdersPerRoot+1:iRoot*nOrdersPerRoot)'));
end
end
end
function x=findzero(n,k,x0,kind)
% Uses Halley's method to find a zero given the starting point x0
% http://en.wikipedia.org/wiki/Halley's_method
ITERATIONS_MAX = 100; % Maximum number of iteration
TOLERANCE_RELATIVE = 1e4; % 16-4 = 12 significant digits
% Setup loop
error = 1;
loopCount = 0;
x = 1; % Initialization value only. It is is not used.
% Begin loop
while abs(error)>eps(x)*TOLERANCE_RELATIVE && loopCount<ITERATIONS_MAX
switch kind
case 1
a = besselj(n,x0);
b = besselj((n+1),x0);
case 2
a = bessely(n,x0);
b = bessely((n+1),x0);
end
xSquared = x0*x0;
error = 2*a*x0*(n*a-b*x0)/(2*b*b*xSquared-a*b*x0*(4*n+1)+(n*(n+1)+xSquared)*a*a);
% Prepare for next loop
x=x0-error;
x0=x;
loopCount=loopCount+1;
end
% Handle maximum iterations
if loopCount>ITERATIONS_MAX-1
warning('Failed to converge to within relative tolerance of %e for n=%f and k=%d in %d iterations', eps(x)*TOLERANCE_RELATIVE, n, k, ITERATIONS_MAX);
end
end
function dJndx = dbesselj(n,x)
% DBESSELJ A function that will generically calculate the
% the derivative of a Bessel function of the first
% kind of order n for all values of x.
%
% Example usage: dJndx = dbesselj(n,x);
%
% INPUT ARGUMENTS
% ================
% n Order of the Bessel function of the first kind
% x Input variable to Bessel function
%
% OUTPUT ARGUMENTS
% ================
% dJndx Derivative of nth order Bessel function of the first
% kind at all values of x
dJndx = n*besselj(n,x)./x - besselj(n+1,x);
end
Javeria
on 18 Aug 2025
Javeria
on 18 Aug 2025
Umar
on 18 Aug 2025
@Javeria,
Yes, your scaling looks correct! Since both your Z-naught-prime and Z-n-prime have units of 1/m, and your eigenvalues chi are properly scaled (roots divided by R), the R-squared over 2 normalization factor is indeed appropriate.
Your numerical verification confirms the analytical formulas are working well (relative errors around 10 to the minus 6 to 10 to the minus 7). The dimensional analysis also checks out: with Z-naught-prime and Z-n-prime having units of inverse meters and your matrix elements being dimensionless as expected, everything appears consistent.
The key point Torsten and I were emphasizing is that you need the "r divided by R" scaling inside all your Bessel function arguments throughout your derivation, which you have correctly implemented in your code where chi equals roots divided by RI.
Thanks to @Torsten for:
*Confirming the integral formulas and catching the sign issue in the J-I integral
*Providing the numerical validation code that helped verify your analytical expressions
*Clarifying that gamma has units of 1/m by definition (roots divided by R)
For next steps:
*Your numerical checks show good agreement (errors around 10^-6), so your formulas are working correctly
*The dimensional analysis looks consistent
*Focus now on making sure your written mathematical derivation reflects the same proper scaling you've implemented in code
This will give you a solid foundation before moving forward with the rest of your analysis.
Javeria
on 18 Aug 2025
You could try like this, but the array E1 (and others) are not defined. So I don't know if it works.
clear all;
close all;
clc;
tic;
tau = [0.4 1];
hold on
for i = 1:numel(tau)
[omega,eta] = compute_eta(tau(i));
plot(2*pi./omega, eta, 'k', 'LineWidth', 1.5); %%% this is T which is T= 2*pi/omega plotting against T
end
hold off
xlabel('$T$', 'Interpreter', 'latex');
ylabel('$|\eta / (iA)|$', 'Interpreter', 'latex');
title('Wave motion amplitude for $R_E = 138$', 'Interpreter', 'latex');
legend({'$\tau=0.4$'}, ...
'Interpreter','latex','Location','northwest');
grid on;
xlim([5 35]); % Match the reference plot
ylim([0 4.5]); % Optional: based on expected peak height
elapsedTime = toc;
disp(['Time consuming = ', num2str(elapsedTime), ' s']);
function [omega,eta0] = compute_eta(tau)
%clear all;
%close all;
%clc;
%tic;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Parameters from your setup
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
N = 100; % Number of N modes
M = 100; % Number of M modes
Q = 100; % Number of Q modes
RI = 34.5; % Inner radius (m)
% ratio = 47.5/34.5; % Ratio R_E / R_I
RE = 47.5; % outer radius (m)
% h = 200/34.5 * RI; % Water depth (m)
h = 200;
d = 38; % Interface depth (m)
b = h-d; % Distance from cylinder bottom to seabed (m)
g = 9.81; % Gravity (m/s^2)
%tau = 1.0; % porosity ratio %%% need loop for this one
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
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% ==============Initialize eta0 array before loop
%===============Initialize omega array before loop
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
eta0 = zeros(size(X_kRE)); % Preallocate
omega = zeros(size(X_kRE)); %%%% omega preallocate
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for idx = 1:length(X_kRE)
wk = X_kRE(idx); % dimensional wavenumber
omega(idx) = sqrt(g * wk * tanh(wk * h)); % dispersion relation
%%%%%%%%%%%%% 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
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,optimset('Display','none'));
%%%%%%%%%%%%% 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
%
% 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
% G(q,1) = -1i * (2*a1(idx)/omega(idx)) ... %%%%% Prof.chen solution
% * ( Z0d * chi(q) * besselj(l, wk*RI) ) ...
% / ( (wk^2 - chi(q)^2) * dbesselj(l, wk*RI) );
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
% G(q,n) = -1i * (2*a1(idx)/omega(idx)) ... %%%%% Prof.chen solution
% * ( Znd(n) * chi(q) * besseli(l, k(n)*RI) ) ...
% / ( (k(n)^2 + chi(q)^2) * dbesseli(l, k(n)*RI) );
end
end
% % % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% %%%%%%%%% 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
% % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% %
% %
%plot(2*pi./omega, 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
% % %
%xlabel('$T$', 'Interpreter', 'latex');
%ylabel('$|\eta / (iA)|$', 'Interpreter', 'latex');
%title('Wave motion amplitude for $R_E = 138$', 'Interpreter', 'latex');
%legend({'$\tau=0.4$'}, ...
% 'Interpreter','latex','Location','northwest');
%grid on;
%xlim([5 35]); % Match the reference plot
%ylim([0 4.5]); % Optional: based on expected peak height
%
%elapsedTime = toc;
%disp(['Time consuming = ', num2str(elapsedTime), ' s']);
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 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 = 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 = 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 = 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
Javeria
on 19 Aug 2025
Categories
Find more on Dates and Time in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
