numerical double integral of an expression with vector defined variables

Question

Nima 2020-4-9

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/516594-numerical-double-integral-of-an-expression-with-vector-defined-variables

评论： darova 2020-4-13

Hi everyone,

I have to compute double integral for the expression G over phi and theta in the following code numerically. Does anyone have any idea how it must be done, since it is not easy to use trapz or integral2 because of the definition of the variables!

c = 4.2;
d = 10;
r3 = linspace(c, d, 100);
theta = linspace(-pi/2, pi/2, 100); 
theta2 = linspace(-pi/2 ,pi/2, 100); 
Phi = linspace(0, 2*pi, 100);
Phi2 = linspace(0, 2*pi, 100);
G =  r3.*(sin(theta2).*cos(theta).*cos(phi2-phi) - cos(theta2).*sin(theta)) ./(r3.^2 + c.^2 - 2.*r3.*c.*(sin(theta2).*sin(theta).*cos(phi2-phi)+cos(theta2).*cos(theta))).^3/2 ;

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

darova 2020-4-9

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/516594-numerical-double-integral-of-an-expression-with-vector-defined-variables#answer_425116

在 MATLAB Online 中打开

Here is a start

phi = linspace(...)             % define value for phi
theta = linspace(...);          % define theta
dphi = phi(2)-phi(1);           % phi step
dtheta = theta(2)-theta(1);     % theta step
[PHI,THETA] = meshgrid(phi,theta);      % create 2D grid
G = PHI .. THETA ...            % calculate G at each point of grid
F = G(1:end-1,1:end-1) + G(1:end-1,2:end) + G(2:end,1:end-1) + G(2:end,2:end);
result = sum(F(:))/4*dphi*dtheta*

Do you know why F variable is so long?

5 个评论
显示 3更早的评论隐藏 3更早的评论

Nima 2020-4-12

在 MATLAB Online 中打开

Many thanks for your answer!

there is a new issue now: I want to compute it for all the points, which means each point at (theta2,phi2) relative to all the points at(theta,phi), for this reason, I have defined the following for loops:

theta = linspace(-pi/2, pi/2, 100);
phi = linspace(0, 2*pi, 100);
r3 = linspace(4.2, 10, 100);
dphi = phi(2)-phi(1);           % phi step
dtheta = theta(2)-theta(1);     % theta step
[PHI,THETA] = meshgrid(phi,theta); 
for i=1:numel(theta)
    for j=1:numel(phi)
        
G{i,j} =   r3.*(sin(THETA).*cos(THETA(i)).*cos(PHI-PHI(j)) - cos(THETA).*sin(THETA(i))) + r3.*sin(THETA).*sin(PHI-PHI(j))./(r3.^2 + c.^2 - 2.*r3.*c.*(sin(THETA).*sin(THETA(i)).*cos(PHI-PHI(j))+cos(THETA).*cos(THETA(i)))).^3/2 ;
    end
end
G{i,j}(isnan(G{i,j})) = 0;
F = G{i,j}(1:end-1,1:end-1) + G{i,j}(1:end-1,2:end) + G{i,j}(2:end,1:end-1) + G{i,j}(2:end,2:end);
result = sum(F(:))/4*dphi*dtheta;

My questions: (1) why result still is a number and not a matrix in despite of G is 100*100 cell?

(2) and must be r3 a vector or matrix?

(3) could you please give a explanation about variable F?

darova 2020-4-12

在 MATLAB Online 中打开

Sorry, big THETA here

F = F .* R^2 .* cos(THETA(2:end,2:end)); % correction of elementary volume

It's late today. Im going to sleep. I'll respond tomorrow

darova 2020-4-13

在 MATLAB Online 中打开

As for volume calculations:

I suggest you to use 3D matrices:

theta = ...
phi = ..
rr = ...
dth = theta(2) - theta(1);
dphi = phi(2) - phi(1);
dr = rr(2) - rr(1);
[T,P,R] = meshgrid(theta,phi,rr);       % create 3D matrices

Integral 3

dV = (your long function) .* R.^2 .*cos(T) *dr*dth*dphi;

Now there is need to 'averaging' as before but for 3D matrices

Maybe there is better idea but i have only this:

dV = dV(1:end-1,1:end-1,1:end-1) + ...
     dV(1:end-1,1:end-1,2:end) + ...
     dV(1:end-1,2:end,1:end-1) + ...
     dV(2:end,2:end,1:end-1) + ...
     % 4 more combinations
dV = sum(dV(:))/8;

There is one more idea, but nor precise one (without averaging). Just create dense mesh

dV1 = dV(2:end,2:end,2:end);
dV1 = sum(dV1);

请先登录，再进行评论。

numerical double integral of an expression with vector defined variables

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

采纳的回答

5 个评论
显示 3更早的评论隐藏 3更早的评论

更多回答（0 个）

另请参阅

类别

标签

产品

Community Treasure Hunt

numerical double integral of an expression with vector defined variables

0 个评论 显示 -2更早的评论隐藏 -2更早的评论

采纳的回答

5 个评论 显示 3更早的评论隐藏 3更早的评论

更多回答（0 个）

另请参阅

类别

标签

产品

Community Treasure Hunt

WeChat

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

5 个评论
显示 3更早的评论隐藏 3更早的评论