用二分法求第一类贝塞尔函数(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.

我们必须使用此代码来查找第一类贝塞尔函数(J0(x))的第n个零。插入一个范围，然后找到我们正在寻找的特定根，这是非常简单的。然而，我们必须绘制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)，.....so on的根。(i=0,1,2，... )

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