Нахождение n-го корня функции Бесселя первого рода (J0 (x)) с помощью деления пополам - PullRequest
Новая
       15

Нахождение n-го корня функции Бесселя первого рода (J0 (x)) с помощью деления пополам

0 голосов
/ 27 сентября 2018

Прежде всего, я бы хотел уточнить, что это задание для школы, поэтому я не ищу решения.Я просто хочу, чтобы меня подтолкнули в правильном направлении.

Теперь о проблеме.

У нас есть код для нахождения корня многочлена с помощью деления пополам: 

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

Я верил, что это сработает, но предоставляемые им корни не в порядке и продолжают повторяться.Проблема, которую я считаю, заключается в моем диапазоне, когда я называю бисекцию2.Я знаю, что [i, j] - не лучший способ сделать это, и я надеялся, что кто-то может направить меня в правильном направлении, как решить эту проблему.

Спасибо.

1 Ответ

0 голосов
/ 27 сентября 2018

Ваша реализация в правильном направлении, но она не совсем верна. 

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].Так как 1-й корень для J0 имеет значение 2,4, корень, который вычисляет ваша функция деления пополам (0,999), на самом деле не является первым корнем.Строки между if ..... end будут проверять, является ли корень, найденный функцией Bisection, на самом деле корнем или нет.Если это так, x_ini примет значение корня, так как следующий корень появится после x = x_ini (или y).

...