利用MATLAB代码实现

好的，以下是MATLAB代码实现：

% 定义被积函数
fun = @(r,theta) besselj(0,k*r)+besselj(1,k*r).*cos(theta);

% 定义积分区间和步长
a = -0.5; b = 0.5;
hx = 0.1; hy = 0.1;

% 计算 x 和 y 方向上的节点数
nx = (b-a)/hx+1; ny = (b-a)/hy+1;

% 生成节点坐标矩阵
[x,y] = meshgrid(linspace(a,b,nx),linspace(a,b,ny));

% 计算每个小矩形的面积
dA = hx*hy;

% 计算每个小矩形的近似积分值（复化梯形公式）
sumf = 0;
for i = 1:nx-1
    for j = 1:ny-1
        sumf = sumf + dA/4*(fun(x(i,j),y(i,j))+fun(x(i,j+1),y(i,j+1))+fun(x(i+1,j),y(i+1,j))+fun(x(i+1,j+1),y(i+1,j+1)));
    end
end
I_trapezoidal = sumf;

% 计算每个小矩形的近似积分值（复化Simpson公式）
sumf = 0;
for i = 1:2:nx-2
    for j = 1:2:ny-2
        sumf = sumf + dA/9*(fun(x(i,j),y(i,j))+4*fun((x(i,j)+x(i+1,j))/2,y(i,j))+fun(x(i+1,j),y(i,j))+4*fun(x(i,j+1),(y(i,j+1)+y(i,j))/2)+2*fun((x(i,j)+x(i+1,j))/2,(y(i,j)+y(i+1,j+1))/2)+4*fun(x(i+1,j+1),(y(i,j+1)+y(i+1,j+1))/2)+fun(x(i,j+1),y(i+1,j+1))+4*fun((x(i,j)+x(i+1,j+1))/2,y(i+1,j+1))+fun(x(i+1,j+1),y(i+1,j+1)));
    end
end
I_Simpson = sumf;

% 输出结果
fprintf('利用复化梯形公式，积分值为：%f\n',I_trapezoidal);
fprintf('利用复化Simpson公式，积分值为：%f\n',I_Simpson);

其中，`k`为常数，可以取10、50或100，`besselj(0,k*r)`和`besselj(1,k*r)`分别表示第一类Bessel函数$J_0$和$J_1$。在实际计算中，我们需要将极坐标形式的被积函数转换成直角坐标形式，这可以通过以下代码实现：

f = @(x,y) fun(sqrt(x.^2+y.^2),atan2(y,x));

然后，在计算每个小矩形的近似积分值时，可以采用上面提到的复化梯形公式和复化Simpson公式。最后，将所有小矩形的近似积分值相加，即可得到整个区域的近似积分值。

希望这个MATLAB代码可以帮助你实现计算。