0

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

现在，对于问题。

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

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] 不是最好的方法，我希望有人能引导我朝着正确的方向解决这个问题。

谢谢你。

4

1 回答 1

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) 实际上并不是第一个根。if.....end 之间的行将检查 Bisection 函数找到的根是否实际上是根。如果是，x_ini 将取根的值，因为下一个根将出现在 x = x_ini（或 y）之后。

于 2018-09-27T07:47:32.630 回答