# How to solve a threedimensional integral of a nested function?

2 views (last 30 days)
Simon on 16 Dec 2013
Commented: Mike Hosea on 7 Jan 2014
Hello,
I'm working with Matlab for a month now at my student research job at university. We are designing a new test rig to optically measure the moistness in a two-phase flow. My task is to implement a calibration curve that associates the relative change of the laser intensity to the measured droplet diameter. The theory behind that method is the Mie-theory or Mie-scattering; just so that you have an idea why I'm doing this.
The function that gives the relative laser intensity for a given droplet diameter is dependent on and needs to be integrated over two angles and a wavelength. The problem is, that I don't know how to implement the integration either with the integral3 command or numerically with a loop.
Using the integral3 command I would need something like
g = @(x,y,z) x.*y.*2.*z;
f = @(x,y,z) 2.*g+3.*y-4.*z;
q = integral3(f,min,max...)
but that is not working and Matlab gives me multiple error messages that are related to the implemented integral3 code.
Another option I thought of was to integrate numerically using a while loop. I splitted every integration variable range into the same number of intervals and calculated the function for every step, hence from f(lambda_0,phi_0,theta_0) to f(lambda_n,phi_n,theta_n). This approach worked pretty well since there also is a sum from 0 to infinity that needs to be calculated throughout the integral.
But in the end when I tried to compute the integral with the trapezium rule it ocurred to me that I do not know the interval width to multipy by since it differs for the integration variables. Or is it generally not allowed to use the trapezium rule for a more-dimensional problem?
Do you have a good idea for me how to solve my problem? If you need I can provide the code I have written so far.
Thank you very much, Simon
Edit: Here is the code I have written so far
J=100;
K=1;
r=3;
R=cos(thetar);
P=zeros(1,J+1); % +1, because Array can't initialize a zeroth position
dP=zeros(1,J+1);
%Initial values for recursive calculation of legendre-polynomials
%BTW is there another chance to calculate the legendre-polynomials in
Matlab? I found the orthpoly function but this is only available in in
the mathematic toolbox
P(0+1)=1;
P(1+1)=R;
P(2+1)=1/2 .* ( 3.*R.^2-1 );
dP(0+1)=0;
dP(1+1)=1;
dP(2+1)=3 .* R;
while r<=J
r=r+1;
P(r)=(((2.*(r-2)+1).*R.*P(r-1)-(r-2).*P(r-2))./((r-2)+1));
dP(r)=((2.*(r-2)+1).*P(r-1)+dP(r-2));
end
s1K=zeros(1,J);
s2K=zeros(1,J);
aK=zeros(1,J);
bK=zeros(1,J);
PSIy = zeros(1,J);
PSIx = zeros(1,J);
CHIx = zeros(1,J);
PSIyI = zeros(1,J);
PSIxI = zeros(1,J);
CHIxI = zeros(1,J);
while K<=J
x = pi*d./lambdar; % d and lambda of the same unit
y = n*x;
PSIy(K)=besselj(y,K);
PSIx(K)=besselj(x,K);
CHIx(K)=PSIx(K)+1i.*bessely(x,K);
PSIyI(K)=(1./2).*(besselj(y,K-1)-besselj(y,K+1));
PSIxI(K)=(1./2).*(besselj(x,K-1)-besselj(x,K+1));
CHIxI(K)=PSIxI(K)+1i.*bessely(x,K);
tanalpha=-(PSIyI.*PSIx-n.*PSIy.*PSIxI)./(PSIyI.*CHIx-n.*PSIy.*CHIxI);
tanbeta=-(n.*PSIyI.*PSIx-PSIy.*PSIxI)./(n.*PSIyI.*CHIx-PSIy.*CHIxI);
aK=tanalpha./(tanalpha-1i);
bK=tanbeta./(tanbeta-1i);
s1K(K)=((2.*K+1)./(K*(K+1))).*(aK(K).*dP(K+1)+bK(K).*(K.*(K+1).*P(K+1)-cos(dP(K+1))));
s2K(K)=((2.*K+1)./(K*(K+1))).*(aK(K).*(K.*(K+1).*P(K+1)-cos(dP(K+1)))+bK(K).*dP(K+1));
K=K+1;
end
f=(((lambdar).^2)/(8*z^2*pi^2)).*(i1.*(sin(phir)).^2+i2.*(cos(phir)).^2);
% This is the function that needs to be integrated by the end of the day
with integration variables lambdar, phir and thetar
i1=abs(s1sum);
i2=abs(s2sum);
s1sum=sum(s1K);
s2sum=sum(s2K);

Mike Hosea on 16 Dec 2013
Edited: Mike Hosea on 16 Dec 2013
The devil here must be in the details. The example you give is easy enough
>> g = @(x,y,z) x.*y.*2.*z;
>> f = @(x,y,z) 2.*g(x,y,z)+3.*y-4.*z;
>> q = integral3(f,0,1,0,1,0,1)
q =
2.775557561562891e-17
Where I have obviously supplied my own xmin, xmax, ymin, ymax, etc. The only important revision here is to correct the definition of f. It is difficult to know whether this was an inadvertent error in a hastily jotted example or indicative of a misunderstanding about how to use function handles. I tend to assume the former, so it would be really helpful to know what problem INTEGRAL3 is having.
Now if your integral has discontinuities internal to the region, then you should be using 'method','iterated'. Your integrand needs to be properly vectorized, i.e. w = f(x,y,z) must return an array w of the same size as x, y, and z (they will be the same size) where w(i,j,k) = f(x(i,j,k),y(i,j,k),z(i,j,k)).
If the integrand has some other internal feedback, then it is not going to work. If you have trouble vectorizing the function because it is too complicated, you can make f work on scalars and then use arrayfun to extend it to array inputs like so
fv = @(x,y,z)arrayfun(f,x,y,z);
integral3(fv,xmin,xmax,ymin,ymax,zmin,zmax);
I would recommend that you really loosen the tolerances up for the first run, a la
integral3(fv,xmin,xmax,ymin,ymax,zmin,zmax,'AbsTol',1e-3,'RelTol',1e-3)
This may be quite slow, unfortunately, since everything is scalarized internally. Note that if you have any singularities, you MUST partition this into multiple integrals, putting the singularities on the boundaries of integrations only. Discontinuities do not have to be on the boundaries, but they really, really (really) should be. If that is a hassle (because the integrand arises from a discretization), then as previously mentioned, you will need to use the iterated method
integral3(fv,xmin,xmax,ymin,ymax,zmin,zmax,'AbsTol',1e-3,'RelTol',1e-3,'method','iterated')
Mike Hosea on 7 Jan 2014
I don't know what i1i2_final does. Presumably d, n, and z are constant scalars that have already been defined (a snapshot of their respective values is taken when each of those three functions is defined.)
Be sure to test with arrays, e.g. f(rand(3,4),rand(3,4),rand(3,4)). If the function only works on scalars, then, yes, you will need arrayfun.

Vaclav Rimal on 16 Dec 2013
Edited: Vaclav Rimal on 16 Dec 2013
First, I suppose ther should be g(x,y,z) instead of g in the definition of function f.
I would solve this by running trapz three times. Before, you need to specify the limits of integration and create corresponding axes, calculate the values of f(x,y,z) for all points using meshgrid.
x = -3:0.01:3;
y = -0.5:0.01:0.5;
z = -1:0.01:1;
[xxx,yyy,zzz]=meshgrid(x,y,z);
i1=trapz(z,f(xxx,yyy,zzz),3); % integration along z
i2=trapz(y,i1); % integration along y
i3=trapz(x,i2); % integration along x
However, creating large 3D array can cause running out of memory. You can avoid meshgrid by including double for loop inside your function, if it helps, like this or something similar:
for i=1:numel(x)
for j=1:numel(y)
f(i,j,:)=func(x(i),y(j),z);
end
end
Simon on 18 Dec 2013
I'm not sure to this point if I can even use a function like trapz to calculate the integral. I'll provide the code I have written so far to give you an idea what it looks like.