Helmholtz problem in circular disk

The following expression

gives the solution for the Helmholtz problem. On the circular disc with center 0 and radius a. For the plot in 3-dimensional graphics of the solutions on Matlab for .

the first row of my data should display the first row of , and ​ as "2.4048, 3.8317, 5.1356

Why do I have this results at (1,1) and (2,1)? WHat should I have to change?

diska = 1;  % Radius of the disk
mmax = 2;   % Maximum value of m
nmax = 2;   % Maximum value of n

% Function to find the k-th zero of the n-th Bessel function
% This function uses a more accurate method for initial guess
besselzero = @(n, k) fzero(@(x) besselj(n, x), [(k-1)*pi, k*pi]);

% Define the eigenvalue k[m, n] based on the zeros of the Bessel function
k = @(m, n) besselzero(n, m);

% Define the functions uc and us using Bessel functions
% These functions represent the radial part of the solution
uc = @(r, t, m, n) cos(n * t) .* besselj(n, k(m, n) * r);
us = @(r, t, m, n) sin(n * t) .* besselj(n, k(m, n) * r);

% Generate data for demonstration
data = zeros(5, 3);
for m = 1:5
    for n = 0:2
        data(m, n+1) = k(m, n);  % Storing the eigenvalues
    end
end

% Display the data
disp(data);

    2.4048         0         0
    5.5201    3.8317    5.1356
    8.6537    7.0156    8.4172
   11.7915   10.1735   11.6198
   14.9309   13.3237   14.7960

% Plotting all in one figure
figure;
plotIndex = 1;
for n = 0:nmax
    for m = 1:mmax
        subplot(nmax + 1, mmax, plotIndex);
        [X, Y] = meshgrid(linspace(-diska, diska, 100), linspace(-diska, diska, 100));
        R = sqrt(X.^2 + Y.^2);
        T = atan2(Y, X);
        Z = uc(R, T, m, n);  % Using uc for plotting
        % Ensure the plot is only within the disk
        Z(R > diska) = NaN;
        mesh(X, Y, Z);
        title(sprintf('uc: n=%d, m=%d', n, m));
        colormap('jet');
        plotIndex = plotIndex + 1;
    end
end

Answer 1

[David Goodmanson]

Hi Athanasios,

the problem is with the table of zeros of the bessel functions. Counting starts with the first nonzero value, so your table of zero locations has an off-by-one issue. I made a quick fix (both your code and the fix do NOT work for n>=3, see below) by replacing

besselzero = @(n, k) fzero(@(x) besselj(n, x), [(k-1)*pi, k*pi]);

with

besselzero = @(n, k) fzero(@(x) besselj(n, x), [(k-(n==0))*pi, (k+1-(n==0))*pi]);

so now the table comes out as

2.4048    3.8317    5.1356
5.5201    7.0156    8.4172
8.6537   10.1735   11.6198
11.7915   13.3237   14.7960
14.9309   16.4706   17.9598

and the plot is

However, as you may be aware, the besselzero function in your code does not work for J_n when n>=3. This is because the search bounds for fzero, i.e. [(k-1)*pi, k*pi], happen to work for n = 0,1,2 but that's it. And use of the fzero function slows things down, although for moderate values of n and m that will not matter much.

An improved function is showed in the code below. It uses Newton's method (which means it can be vectorized) and works for any reasonable n and m. Calculating e.g. the first 100 roots of J_6 takes about one millisecond.

function x = bessel0j(n,q,opt)

% a row vector of the first q roots of bessel function Jn(x), integer order.
% if opt = 'd', row vector of the first q roots of dJn(x)/dx, integer order. 
% if opt is not provided, the default is zeros of Jn(x).
% all roots are positive, except when n=0,
% x=0 is included as a root of dJ0(x)/dx (standard convention).
%
% starting point for for zeros of Jn was borrowed from Cleve Moler, 
% but the starting points for both Jn and Jn' can be found in
% Abramowitz and Stegun 9.5.12, 9.5.13. 
%
% function x = bessel0j(n,q,opt)

% David Goodmanson

k = 1:q;
if nargin==3 & opt=='d'
    beta = (k + n/2 - 3/4)*pi;
    mu = 4*n^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(n,x)./ ...
            (besselj(n,x).*((n^2./x.^2)-1) -besseljd(n,x)./x);
        x = xnew;    
    end
    if n==0
        x(1) = 0;    % correct a small numerical difference from 0
    end
else
    beta = (k + n/2 - 1/4)*pi;
    mu = 4*n^2;
    x = beta - (mu-1)./(8*beta) - 4*(mu-1)*(7*mu-31)./(3*(8*beta).^3);
    for j=1:8
        xnew = x - besselj(n,x)./besseljd(n,x);
        x = xnew;
    end
end

[Athanasios Paraskevopoulos]

Thank you very much for your help

Answer 2

[Torsten]

According to your code, data(i,j) is a zero of besselj(j-1,x) in the interval [(i-1)*pi i*pi].

So let's plot besselj(1,x) and besselj(2,x) in the interval [0 pi]:

x = linspace(0,pi,100);
hold on
plot(x,besselj(1,x))
plot(x,besselj(2,x))
hold off
grid on

Thus data(1,2) = data(1,3) = 0 is correct.