使用平分法找到第一种贝塞尔函数 (J0(x)) 的第 n 个根

Question

首先，我只想澄清一下，这是学校的作业，所以我不是在寻找解决方案。我只是想被推向正确的方向。

现在，对于问题。

我们有使用二分法查找多项式根的代码：

function [root, niter, rlist] = bisection2( func, xint, tol )
% BISECTION2: Bisection algorithm for solving a nonlinear equation
%             (non-recursive).
%
% Sample usage:
%   [root, niter, rlist] = bisection2( func, xint, tol )
% 
%  Input:
%     func    - function to be solved
%     xint    - interval [xleft,xright] bracketing the root
%     tol     - convergence tolerance (OPTIONAL, defaults to 1e-6)
%
%  Output:
%     root    - final estimate of the root
%     niter   - number of iterations needed  
%     rlist   - list of midpoint values obtained in each iteration.
% First, do some error checking on parameters.
if nargin < 2
fprintf( 1, 'BISECTION2: must be called with at least two arguments' );
error( 'Usage:  [root, niter, rlist] = bisection( func, xint, [tol])' );
end
if length(xint) ~= 2, error( 'Parameter ''xint'' must be a vector of length 2.' ), end  
if nargin < 3, tol = 1e-6; end
% fcnchk(...) allows a string function to be sent as a parameter, and
% coverts it to the correct type to allow evaluation by feval().
func = fcnchk( func );
done  = 0;
rlist = [xint(1); xint(2)];
niter = 0;
while ~done
% The next line is a more accurate way of computing
% xmid = (x(1) + x(2)) / 2 that avoids cancellation error.
xmid = xint(1) + (xint(2) - xint(1)) / 2;
fmid = feval(func,xmid);
if fmid * feval(func,xint(1)) < 0
xint(2) = xmid;
else
xint(1) = xmid;
end
rlist = [rlist; xmid];
niter = niter + 1;
if abs(xint(2)-xint(1)) < 2*tol || abs(fmid) < tol
done = 1;
end
end
root = xmid;
%END bisection2.

我们必须使用此代码来查找第一种贝塞尔函数的第 n 个零 (J0(x((。插入一个范围，然后找到我们正在寻找的特定根非常简单。但是，我们必须绘制 Xn 与 n 的关系图，为此，我们需要能够计算出与 n 相关的大量根。为此，我编写了以下代码：

bound = 1000;
x = linspace(0, bound, 1000);
for i=0:bound
for j=1:bound
y = bisection2(@(x) besselj(0,x), [i,j], 1e-6)
end
end

我相信这会奏效，但它提供的根源并不有序，并且不断重复。我认为问题是我调用 bisection2 时的范围。我知道[i，j]不是最好的方法，并希望有人能引导我朝着如何解决这个问题的正确方向前进。

谢谢。

Answer 1

你的实现方向是正确的，但它并不完全正确。

bound = 1000;
% x = linspace(0, bound, 1000);  No need of this line.
x_ini = 0; n =1; 
Root = zeros(bound+1,100);  % Assuming there are 100 roots in [1,1000] range
for i=0:bound                                            
for j=1:bound                                          
y = bisection2(@(x) besselj(i,x), [x_ini,j], 1e-6);   % besselj(i,x) = Ji(x); i = 0,1,2,3,...
if abs(bessel(i,y)) < 1e-6
x_ini = y;                                         % Finds nth root
Root(i+1,n) = y;
n = n+1;
end
end
n = 1;
end

我已经将代码中的 besselj(0，x( 替换为 besselj(i，x(。这不仅为您提供了J0(x(的根，而且还为您提供了J1(x(，J2(x(，J3(x(的根，.....等等。(i=0,1,2,... (

我在您的代码中所做的另一个更改是将 [i，j] 替换为 [x_ini，j]。最初x_ini=0 & j=1。这将尝试在区间 [0,1] 中查找根。由于 J0 的第一个根出现在 2.4 处，因此平分函数计算的根 (0.999( 实际上并不是第一个根。如果.....end将检查平分函数找到的根是否实际上是根。如果是，x_ini将获取根的值，因为下一个根将在 x = x_ini(或 y(之后出现。