کد متلب روش های عددی، مثال های کتاب سادیکو، فصل 3، مثال 3-10

clear all

% Compare MATALB Bessel function to a customly created Bessel function

% Inputs:

x = 0.1:0.1:2; % input to bessel function

m = 0; %order of bessel function

n = 6;      % Newton-Cotes rule

N = 840;  % Newton-Cotes N  

% Compare with built in MATLAB

J = besselj(m,x);

% Create custom Bessel function

lwr = 0;    %lower limit of integration

uppr= pi;   %upper limit of integration

theta = linspace(0,pi,N+1); 

for ii = 1:length(x)    

    % Create function to be integrated

    A = cos(x(ii)*sin(theta)-m*theta)/pi;    

    %Integrate with newton-cotes method

    J_m(ii) = nc_method(A,lwr,uppr,n);

end

% Difference between MATLAB and custom Bessel

compare = J-J_m;

% eps is the Spacing of floating point numbers in MATLAB.

% Any number smaller than in magnitude than eps is essentially zero

% Try 1 + eps, this is equal to 1.

v_eps = [eps, eps]; 

x_eps = [x(1),x(end)];

figure(1),

    subplot(211),plot(x,J,x,J_m,':')

        title('Comparison between MATLAB and Custom Bessel function')

        legend('MATLAB','Custom')

    subplot(212),plot(x,compare,'-o',x_eps,v_eps,x_eps,-v_eps)

        title('Difference between MATLAB and Custom Bessel')

        legend('difference','+eps','-eps')

%figure(2),plot(theta,A)