numerical double integral of an expression with vector defined variables

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 ;

Answer 1

[darova]

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?

[Nima]

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]

I don't understand your new issue. Can you explain? Are you trying to calculate 3 integral now?

I assumed that you are trying to find volume under the surface. I simply replaced Z value with mean value of 4 neighbour vertices. So elementary volume is

I didn't notice before that you are trying to calculate something in spherical system. Elementary volume is spherical as following (volume at the top ( ) is small and volume at is bigger)

So my previous code should be like this

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);
F = F .* R^2 .* cos(theta(2:end,2:end));    % correction of elementary volume
result = sum(F(:))/4*dphi*dtheta*

[Nima]

yes, you got the point correctly, but is this true to do so:

result = sum(F(:))/4*dphi*dtheta.*r.^2.*sin(theta);

because, when do like

F = F .* R^2 .* cos(theta(2:end,2:end)); 

then

Error using  .* 
Matrix dimensions must agree.

and also based on the description, cos or sin in your cod?

[darova]

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]

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);