Clear Filters
Clear Filters

Helmholtz problem in circular disk

9 views (last 30 days)
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
% Display the 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
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));
plotIndex = plotIndex + 1;

Accepted Answer

David Goodmanson
David Goodmanson on 16 Mar 2024
Edited: David Goodmanson on 17 Mar 2024
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]);
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;
if n==0
x(1) = 0; % correct a small numerical difference from 0
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;

More Answers (1)

Torsten on 16 Mar 2024
Moved: Torsten on 16 Mar 2024
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
hold off
grid on
Thus data(1,2) = data(1,3) = 0 is correct.

Community Treasure Hunt

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

Start Hunting!