对一个函数进行三重积分，函数在附件中。

在球坐标系中对这个式子积分，按照剖分定义不断叠加for循环，因剖分次数太多，所以电脑一直出不来结果。想利用矩阵化代替for，写了以后出现了问题。贴出原程序和矩阵化以后的程序，想请大家帮我看看下一步要怎么进行。clear

clc

bx=0;by=0;

b_x=0.05;b_y=0.05;

Nx=49;Ny=49;

f1=1200;

lamda1=340/f1;

k1=2*pi/lamda1;

f2=1e9;

lamda2=3e8/f2;

k2=2*pi/lamda2;

Pt=1000;%w

Pt0=10*log(Pt);

thetav=0.75*pi/180;%天线波束宽度

ge=40;%dB,增益

A0=1;

Z0=A0^2/(4*pi*Pt*ge);%阻抗

%%%%%%%%%%%%%%%环境因子

T=288.15;  %K

e=10.02;   %hpa

P=1013.25;  %hpa

W=10000;      % w   功率

rho_0=1.29; % Kg/m^3   空气密度

c_0=331.6+0.6*15.15;    % 声音速度  温度15.15摄氏度

r_0=5*10^(-3);          % 单个声源半径

A=5*10^(-4);            % 声波传播路径上的衰减

%%%%%%%%%%%%%%%%

theta2=linspace(2*pi/180,179*pi/180,300);%散射角

% theta2=[179*pi/180,119*pi/180,59*pi/180,9*pi/180];

H0=10;%不均匀体内参考点高度

distance=2*H0./tan(theta2/2);%距离

R0=(distance/2)./cos(theta2/2);%波传播“标准路径”

gama=(distance/2)/(2*R0);

radius=R0*thetav/(4*sin(2*gama));%球域的半径

H=H0+radius;%球中心与坐标原点的距离

Z=radius./H;

fei1=linspace(0,2*pi,200);

%%%%%%%%%%%%%%%%

h=2*pi/200;

F=zeros(1,300);

E01=zeros(1,300);

Q1=[];

Q2=[];

Q3=[];

Q4=[];

G=[];

for m=1:length(theta2)

theta1=linspace(0,asin(Z(m)),200);

g=asin(Z(m))/200;

for l=1:length(theta1)

V1=H(m)*cos(theta1(l))-sqrt(radius(m)^2-H(m)^2*(sin(theta1(l))^2));

V2=H(m)*cos(theta1(l))+sqrt(radius(m)^2-H(m)^2*(sin(theta1(l))^2));

r=linspace(V1,V2,200);

t=(V2-V1)/200;

for v=1:length(fei1)

for k=1:200

fei0=v*h;

theta0=l*g;

r0=V1+k*t;

distance0=distance(m)/2;%发收点与坐标系原点的距离

if  fei0>0&&fei0

R1=distance0^2+r0^2-2*distance0*r0*cos(pi/2+theta2(m)/2);

R2=distance0^2+r0^2-2*distance0*r0*cos(pi/2-theta2(m)/2);

else

R1=distance0^2+r0^2-2*distance0*r0*cos(pi/2-theta2(m)/2);

R2=distance0^2+r0^2-2*distance0*r0*cos(theta2(m)/2+pi/2);

end

%%%函数

if theta0==0

theta0=eps;

if fei0==pi/2||fei0==3*pi/2

fei0=fei0+eps;

end

end

S_P=sin(fei0).*sin(Nx*pi*b_x.*sin(theta0).*cos(fei0)/lamda1).*sin(Ny*pi*b_y.*sin(theta0).*cos(fei0)/lamda1)./(sin(pi*b_x.*sin(theta0).*cos(fei0)/lamda1).*sin(pi*b_y.*sin(theta0).*cos(fei0)/lamda1));

U=(r0*cos(fei0)-H0)/(sqrt(r0.^2+H0.^2-2*r0*H0*cos(theta0)));%夹角余弦值

Q1=exp(-1i*2*k2.*sin(theta2(m)/2)*U.*sqrt((H(m)-radius(m)).^2+r0^2-2.*(H(m)-radius(m)).*r0.*cos(theta0)));

for mm=1:Nx

for nn=1:Ny

x=-(Nx+1)/2*b_x+mm*b_x;

y=-(Ny+1)/2*b_y+nn*b_y;

Q2=exp(1i*2*k2.*log(r0)*(x*cos(theta0).*cos(fei0)+y*cos(theta0).*sin(fei0)));

end

end

Q3=exp(-A*r0)*r0^2;

Q4=1/(R1*R2);

G=S_P*Q1*Q2*Q3*Q4;

F(m)=F(m)+G*g*h*t;

Q5=exp(-1i*k2*R1)/R1;

E01(m)=E01(m)+Q5;

end

end

end

end

distance0=distance/2;

Es0=(0.184*10^-6)*k2^2*sqrt(W*rho_0*c_0/(2*pi*r_0^2)).*exp(-1i*k2*2*sqrt(H0.^2+distance0.^2))./(4*pi).*F;

Es=sqrt(30*Pt*ge).*Es0;

E0=sqrt(30*Pt*ge).*E01;

XX=abs(E0);

YY=abs(Es);

ZZ=angle(Es);

下面是我修改后的

clear

clc

bx=0;by=0;

b_x=0.05;b_y=0.05;

Nx=49;Ny=49;

f1=1200;

lamda1=340/f1;

k1=2*pi/lamda1;

f2=1e9;

lamda2=3e8/f2;

k2=2*pi/lamda2;

Pt=1000;%w

Pt0=10*log(Pt);

thetav=0.75*pi/180;%天线波束宽度

ge=40;%dB,增益

A0=1;

Z0=A0^2/(4*pi*Pt*ge);%阻抗

%%%%%%%%%%%%%%%环境因子

T=288.15;  %K

e=10.02;   %hpa

P=1013.25;  %hpa

W=10000;      % w

rho_0=1.29; % Kg/m^3   空气密度

c_0=331.6+0.6*15.15;    % 声音速度  温度15.15摄氏度

r_0=5*10^(-3);          % 单个声源半径

A=5*10^(-4);            % 声波传播路径上的衰减

%%%%%%%%%%%%%%%%

theta2=linspace(2*pi/180,179*pi/180,10);%散射角

% theta2=[179*pi/180,119*pi/180,59*pi/180,9*pi/180];

H0=10;%不均匀体内参考点高度

distance=2*H0./tan(theta2/2);%距离

R0=(distance/2)./cos(theta2/2);%波传播“标准路径”

gama=(distance/2)/(2*R0);

radius=R0*thetav/(4*sin(2*gama));%球域的半径

% V0=4*pi*radius.^3/3;%散射体体积

% length=lamda2./(2*sin(theta2/2));%不均匀体尺寸

H=H0+radius;%球中心与坐标原点的距离

Z=radius./H;

fei1=linspace(0,2*pi,10);

%%%%%%%%%%%%%%%%

h=2*pi/10;

F=zeros(1,10);

E01=zeros(1,10);

Q1=[];

Q2=[];

Q3=[];

Q4=[];

G=[];

A=zeros(1,10);

B=[];

C=[];

D=[];

E=[];

for m=1:length(theta2)

theta1=asin(Z(m));

B=[B,theta1];

for l=1:length(theta1)

V1=H(m)*cos(theta1(l))-sqrt(radius(m)^2-H(m)^2*(sin(theta1(l))^2));

V2=H(m)*cos(theta1(l))+sqrt(radius(m)^2-H(m)^2*(sin(theta1(l))^2));

C=[C,V1(l)];

D=[D,V2(l)];

E=[A B C D];

%             for k=1:200

%                 fei0=v*h;

%                 theta0=l*g;

%                 r0=V1+k*t;

%                 distance0=distance(m)/2;%发收点与坐标系原点的距离

%                 if  fei0>0&&fei0

%                 R1=distance0^2+r0^2-2*distance0*r0*cos(pi/2+theta2(m)/2);%发射波传播实时路径

%                 R2=distance0^2+r0^2-2*distance0*r0*cos(pi/2-theta2(m)/2);%接收波传播实时路径

%                 else

%                     R1=distance0^2+r0^2-2*distance0*r0*cos(pi/2-theta2(m)/2);%发射波传播实时路径

%                     R2=distance0^2+r0^2-2*distance0*r0*cos(theta2(m)/2+pi/2);%接收波传播实时路径

%                 end

%                 %%%函数

%                 if theta0==0

%                         theta0=eps;

%                         if fei0==pi/2||fei0==3*pi/2

%                            fei0=fei0+eps;

%                         end

%                 end

S_P=@(theta,fei) sin(fei).*sin(Nx*pi*b_x.*sin(theta).*cos(fei)/lamda1).*sin(Ny*pi*b_y.*sin(theta).*cos(fei)/lamda1)./(sin(pi*b_x.*sin(theta).*cos(fei)/lamda1).*sin(pi*b_y.*sin(theta).*cos(fei)/lamda1));

U=@(theta,fei,r) (r*cos(fei)-H0)/(sqrt(r.^2+H0.^2-2*r*H0*cos(theta)));%夹角余弦值

Q1=@(theta,r)exp(-1i*2*k2.*sin(theta2(m)/2)*U.*sqrt((H(m)-radius(m)).^2+r^2-2.*(H(m)-radius(m)).*r.*cos(theta)));

for mm=1:Nx

for nn=1:Ny

x=-(Nx+1)/2*b_x+mm*b_x;

y=-(Ny+1)/2*b_y+nn*b_y;

Q2=@(theta,fei,r) exp(1i*2*k2.*log(r)*(x*cos(theta).*cos(fei)+y*cos(theta).*sin(fei)));

end

end

Q3=@(r) exp(-A*r)*r^2;

Q4=1/(R1*R2);

G=S_P*Q1*Q2*Q3*Q4;

a=A(1,m);

b=B(1,m);

c=C(1,m);

d=D(1,m);

F=dblquad(G,a,b,c,d,1.0e-3);

for v=1:length(fei1)

F=F*g;

end

%                 F(m)=F(m)+G*g*h*t;

Q5=exp(-1i*k2*R1)/R1;

E01(m)=E01(m)+Q5;

end

end

我的第一个问题是，fei积分上下限是定值，r和theta则是变量，所以我把三重积分企图化成二重再在外面一重，可是我的函数本身是r,theta,fei的函数，这样拆开积分可以进行吗？如果可以要怎么拆呢？

第二个问题是，被积函数的分母R1，R2是关于实时的r0的函数，并且有fei0所导向。在没有求实时的 fei0、theta0、r0的情况下R1和R2需要怎么处理？

我接触matlab时间不久，有很多方面都还不足，可能这些问题都很简单，但还是希望大家能多多指出我的问题。

2017-12-22 22:49 上传

点击文件名下载附件

16.58 KB, 下载次数: 4