Hi all.
These questions may be the result of mathematical ignorance regarding Bessel
functions, but anyway...
Is it possible to evaluate Bessel functions of complex argument in Maxima?
It appears not for bessel_j as shown here, but I thought I might check with
the members of this list:
(%i68) bessel_j(0,1.0+0.01*%i);
(%o68) BESSEL_J(0, 0.01 %I + 1.0)
(%i69) bessel_j(0,1.0);
(%o69) 0.76519768655797
Bessel_k on the other hand, seems to be evaluable
(%i72) bessel_k(0,1.0);
(%o72) 0.42102443824071
(%i73) bessel_k(0,1.0+0.1*%i);
(%o73) 0.41593669045619 - 0.05982168051612 %I
I would like to examine the behaviour of a complicated expression (deltan1)
arising in the following Maxima batch file, over a range of frequencies
(freq) and wavenumbers (k). I'm interested in places where the expression
goes towards infinity as values of the variables freq and k are altered.
As soon as the Bessel function arguments go complex my scheme goes belly up,
so suggestions about how to go about getting numerical and graphical output
pertinent to that examination with Maxima are also welcome. Alternatively,
how to simplify the expression to a point where my mind's eye can take it
all in?
Putting aside that issue, I'm also wondering whether anyone has comments on
any pitfalls regarding the numerical capabilites of Maxima I may have
overlooked in this batch file.
With thanks for your interest,
Mike Thomas.
/* Parameters */
a: 0.2$ /* Borehole diameter, metres */
alpha_1: 1500$ /* Sound velocity in liquid, metres/sec*/
rho_1: 1000$ /* Fluid density, kg/m^3*/
alpha_2: 1830$ /* Compressional velocity in solid 1, meters/sec */
beta_2: 914$ /* Shear velocity in solid 1, metres/sec */
rho_2: 2000$ /* Density of solid 1, kg/m^3 */
/* freq: 2500$ */ /* Excitation frequency, Hz */
/* k: 300$ */ /* Wavenumber */
omega: 2 * %pi * freq$
v1: (k^2 - omega^2/(alpha_1^2))^(1/2)$
v2: (k^2 - omega^2/(alpha_2^2))^(1/2)$
v3: (k^2 - omega^2/(beta_2^2))^(1/2)$
h: float ( rho_1/rho_2 * ( a * omega / beta_2 )^2 )$
det111: -2 * a * v2 * bessel_k( n+1, a*v2 ) - bessel_k ( n, a * v2 ) * (2 *
n * (n-1) + a * a * (2*k*k - omega*omega / (beta_2*beta_2) ))$
det112: 2 * n * (-a * v3 * bessel_k ( n+1, a * v3 ) + (n-1) * bessel_k ( n,
a * v3))$
det113: 2 * k * a * (a *v3 * bessel_k ( n+1, a * v3 ) + bessel_k ( n, a *
v3 ) * (n * (n-1) + (a * v3)^2))$
det121: 2 * n * (-a * v2 * bessel_k ( n+1, a * v2 ) + (n-1) * bessel_k ( n,
a * v2))$
det122: -2 * a * v3 * bessel_k ( n+1, a * v3 ) - bessel_k ( n, a * v3 ) * (2
* n * (n-1) + (a * v3)^2)$
det123: 2 * a * k * n * (a * v3 * bessel_k ( n+1, a * v3 ) - (n-1) *
bessel_k ( n, a * v3 ))$
det131: 2 * k * a * (-a * v2 * bessel_k ( n+1, a * v2 ) + n * bessel_k ( n,
a * v2 ))$
det132: -a * k * n * bessel_k ( n, a * v3 )$
det133: a * a * (omega*omega / (beta_2*beta_2) - 2 * k * k) * ( - a * v3 *
bessel_k ( n+1, a * v3 ) + n * bessel_k ( n, a * v3 ) )$
det211: a * v2 * bessel_k ( n+1, a * v2 ) - n * bessel_k ( n, a * v2 )$
det212: n * bessel_k ( n, a * v3 )$
det213: k * a * (-a * v3 * bessel_k ( n+1, a * v3 ) + n * bessel_k ( n, a *
v3 ))$
det221: 2 * n * (-a * v2 * bessel_k ( n+1, a * v2 ) + (n-1) * bessel_k ( n,
a * v2 ))$
det222: -2 * a * v3 * bessel_k ( n+1, a * v3 ) - bessel_k ( n, a * v3 ) * (2
* n * (n-1) + (a * v3)^2)$
det223: 2 * a * k * n * (a * v3 * bessel_k ( n+1, a * v3 ) - (n-1) *
bessel_k ( n, a * v3 ))$
det231: 2 * k * a * (-a * v2 * bessel_k ( n+1, a * v2 ) + n * bessel_k ( n,
a * v2 ))$
det232: -a * k * n * bessel_k ( n, a * v3 )$
det233: a * a * (omega*omega / (beta_2*beta_2) - 2 * k * k) * ( a * v3 *
bessel_k ( n+1, a * v3 ) + n * bessel_k ( n, a * v3 ) )$
mat1: matrix ( [det111, det112, det113],
[det121, det122, det123],
[det131, det132, det133] )$
det1: determinant ( mat1 )$
mat2: matrix ( [det211, det212, det213],
[det221, det222, det223],
[det231, det232, det233] )$
det2: determinant ( mat2 )$
deltan: det1 * ( v1 * a * bessel_j ( n+1, v1 * a ) + n * bessel_j ( n, v1
* a ) ) + h * det2 * bessel_j ( n, v1 * a )$
deltan1: deltan,n=1;
plot3d(deltan1,[freq,100,2000],[k,15,50]);
/* s1: solve ( deltan, k ), n=1$ */