How to create a trapezoidal rule routine to evaluate the following Bessel integral without Matlab functions. ME 1101

Answer 1

[Christopher Salerno]

clear all
close all
clc

% Integrates & plots Bessel functions of the first kind
a = 0; % x intial
b = 20;  % x final
c = 0;  % theta intial
d = pi;  % theta final
n = 1;  % order of the Bessel integral
int = 100;  % intervals for x and theta

theta = linspace(c,d,101);
x = linspace(a,b,101);

for i = 1:length(x)

    J_x = (1/pi)*cos(n*theta-x(i)*sin(theta));
    J(i) = (d-c)*(sum(J_x)+sum(J_x(2:length(J_x)-1)))/(2*int);
   
end

plot(x,J)
grid on
grid minor
title('Bessel Function of the First Kind')
xlabel('x')
ylabel('Jn(x)')

% Check with the Bessel function
y = besselj(1,x);
figure
plot(x,y,'b')
grid on
grid minor
title('Bessel Function of the First Kind')
xlabel('x')
ylabel('J_{1}(x)')