【求助】怎么用matlab画带有贝塞尔函数的三维积分图

本帖最后由 z晓晓 于 2021-10-13 21:16 编辑

麻烦请问一下，有大佬知道我的程序错在哪里吗？我想利用的是图2的公式（2），画图3的第一张图。但是matlab一直报错，如图一。我改了很多次了仍然解决不了。希望有朋友可以帮忙解答.>> lamda=6328e-10;
>>R=0.01;
>>n=1.5;
>>fai=pi/180;
>>k=2*pi/lamda;
>>i=sqrt(-1);
>>  r1=linspace(-0.5,0.5);
>>z=linspace(0,1000);
>>[R1,Z]=meshgrid(r1,z);
>>  E2=@(r) besselj(0,k*r*R1./Z)*exp(i*k*(r^2+R1.^2)/2/Z)*exp(-i*k*(n-1)*fai.*r)*r;
>>E3=integral(E2,0,R);

本帖最后由 zzz在路上 于 2021-10-14 08:13 编辑

报错原因是r和z的维数不一致导致的错误。这个问题的思路应该是用for循环：
>>z=linspace(0,1000);
>>[R1,Z]=meshgrid(r1,z);
>>  E2=@(r) besselj(0,k*r*R1./Z)*exp(i*k*(r^2+R1.^2)/2/Z)*exp(-i*k*(n-1)*fai.*r)*r;
>>E3=integral(E2,0,R);
改为：
for index_z=1:length(z)           
      for index_r=1:length(r1)
           E2(index_r,index_z)=@(r) besselj(0,k*r*r1(index_r)./z(index_z))*exp(i*k*(r^2+r1(index_r).^2)/2/z(index_z))*exp(-i*k*(n-1)*fai.*r)*r;
           E3(index_r,index_z)=integral(E2,0,R);
     end
end



未详细验证，大体思路是这个

本帖最后由 一模 于 2021-10-23 23:16 编辑

我感觉是因为函数表达式里的r是一个值，而程序里的R1,Z都是矩阵，所以会造成矩阵纬度不匹配的问题，楼上的大佬说得对，可以用for循环解决，我改了一下程序。
clear
lamda=6328e-10;
R=0.01;
n=1.5;
fai=pi/180;
k=2*pi/lamda;
i=sqrt(-1);
r1=linspace(-0.5,0.5);
z=linspace(0,1000);
% [R1,Z]=meshgrid(r1,z);
%   E2=@(r) besselj(0,k.*r.*R1./Z).*exp(i*k.*(r^2+R1.^2)/2/Z).*exp(-i*k.*(n-1).*fai.*r).*r;
% E2=@(r) besselj(0,k.*r.*0.5./1).*exp(i*k*(r.^2+0.5^2)/2/1);
% E3=integral(E2,0,R);
E3=zeros(length(r1),length(z));
for m=1:length(r1)
    for n=1:length(z)
%     E2=@(r) besselj(0,k.*r.*r1(m)./z(n)).*exp(i*k*(r.^2+r1(m).^2)/2/z(n))
    E2=@(r) besselj(0,k.*r.*r1(m)./z(n)).*exp(i*k*(r.^2+r1(m).^2)/2/z(n)).*exp(-i*k.*(n-1).*fai.*r).*r;
    E3(m,n)=abs(integral(E2,0,R));
    end
end
mesh(E3)
因为我的电脑太慢了，所以直接假设透过率函数为1进行了测试，跑通了，不知道是不是题主想要的。for循环的效率比较低，如果题主有更高效率的方法，我们可以继续交流一下