How to solve a threedimensional integral of a nested function?

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

Answer 1

[Mike Hosea]

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

[Simon]

Well, I just wrote g and not g(x,y,z) in my definition of f, that's why the error occured. So now I' am also able to solve the integral like you did. But I am still not sure if I can solve the integral like this. I think it would be best if I provide (parts of) the code I have written so far so that you can see what I'm struggling with. But I really thank you so far!

Edit: When I'm trying to do this:

[ i1,i2 ] = i1i2_final( thetar,lambdar,d,n );

f = @(lambdar,phir,thetar) (((lambdar).^2)./(8.*z.^2.*pi.^2)).*(i1( thetar,lambdar,d,n ).*(sin(phir)).^2+i2( thetar,lambdar,d,n ).*(cos(phir)).^2);

q = integral3(f,lambda_0,lambda_h,phi_0,phi_h,theta_0,theta_h);

I get the following error messages:

Subscript indices must either be real positive integers or logicals.

Error in @(lambdar,phir,thetar)(((lambdar).^2)./(8.*z.^2.*pi.^2)).*(i1(thetar,lambdar,d,n).*(sin(phir)).^2+i2(thetar,lambdar,d,n).*(cos(phir)).^2)

Error in integral3>@(y,z)FUN(x(1)*ones(size(z)),y,z) (line 139) @(y,z)FUN(x(1)*ones(size(z)),y,z), ...

Error in integral2Calc>integral2t/tensor (line 229) Z = FUN(X,Y); NFE = NFE + 1;

Error in integral2Calc>integral2t (line 56) [Qsub,esub] = tensor(thetaL,thetaR,phiB,phiT);

Error in integral2Calc (line 10) [q,errbnd] = integral2t(fun,xmin,xmax,ymin,ymax,optionstruct);

Error in integral3/innerintegral (line 138) Q1 = integral2Calc( ...

Error in integralCalc/iterateScalarValued (line 314) fx = FUN(t);

Error in integralCalc/vadapt (line 133) [q,errbnd] = iterateScalarValued(u,tinterval,pathlen);

Error in integralCalc (line 76) [q,errbnd] = vadapt(@AtoBInvTransform,interval);

Error in integral3 (line 122) Q = integralCalc(@innerintegral,xmin,xmax,integralOptions);

[Mike Hosea]

Seems like we're jumping into the middle here. Let's focus on getting the integrand function f defined in MATLAB such that for any x, y, and z, you can evaluate f(x,y,z) and get the right number back. We aren't ready to think about integrating anything until we can define the integrand, plug in values for the inputs, and observe that we are getting the correct output for those inputs. The error message suggests to me that you are thinking of i1 as a function, but you have instead defined it as value or array of values. MATLAB is not a symbolic environment. It doesn't maintain the relationships of undefined variables as you string expressions along. If you say z = x + y, then x and y had better already be defined as values, and MATLAB will add them and store the result in z. The variable z then has no connection whatsoever to x and y. It is not a function of x and y, just some value or array of values that was once computed with then-current values of those variables. To make a connection, you must create an anonymous function. Instead of z = x + y, you would write z = @(x,y)x + y, where the prefix @(x,y) just means "The expression that follows is a function of x and y." Then to evaluate z, you supply x and y, e.g. z(1,pi) = 4.14159... You can also define z = @(x,y)x + y + c, and then c had better be defined at that moment. The resulting function z(x,y) will use the snapshot of c that was available at the moment z was defined. Once you sort out how to define your integrand f so that you can pass in scalars for x, y, and z, then you can vectorize it for INTEGRAL3 using

fv = @(x,y,z)arrayfun(f,x,y,z)

ARRAYFUN will make f work when x, y, and z are arrays, evaluating f(x(i),y(i),z(i)) to produce each element of fv(x,y,z). Once we have fv in hand, then, and only then, will we be able to move on to any integration issues.

[Simon]

First of all, I hope you had a nice holiday and a good start into 2014. I was on vacation myself and had no chance to work on the problem.

What you mentioned above could be the reason for all the error messages I recieve. I tried to implement a function for i1 and i2 so that I can define f as:

    [ i1 ] = @(theta,lambda) i1i2_final( theta,lambda,d,n );

    [ i2 ] = @(theta,lambda) i1i2_final( theta,lambda,d,n );

    f = @(lambda,phi,theta) (((lambda).^2)./(8.*z.^2.*pi.^2)).*(i1(theta,lambda).*(sin(phi)).^2+i2(theta,lambda).*(cos(phi)).^2);

The definitions themselves seem to work but when I'm trying to integrate " f " I get some error messages again. Any ideas what I did wrong this time?

Do I have to use ARRAYFUN at this point to finally integrate the function?

And I have defined f the way, that I can calculate a value for any x,y,z so that I could create a three-dimensional array with function values ranging from xmin,ymin,zmin to xmax,ymax,zmax. But even if I would do that, how could I integrate over this array?

Thank you for your help!

[Mike Hosea]

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.

Answer 2

[Vaclav Rimal]

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]

The problem is that my function is really complicated.

To give you an idea, it looks somewhat like this:

I / I_0 (x,y,z) = G(x,y) * (x*constant* (s1(x,y)*f(z)+s2(x,y)*f(z)) )

And the functions s1 and s2 are a sum theirselves and a recurrence formula in them. So I'm not sure if it will work this way.

Any idea on how to implement a Monte Carlo integration?

[Vaclav Rimal]

In my opinion, in order to reliably estimate the integral you have to evaluate the function in many points in the xyz space anyway, independently on the integration method you choose. So, how complicated your function is shouldn't matter and I can't see how M-C integration would avoid it. Or is this the problem, that you cannot calculate the function in lots of points in a finite time?

What can be an issues is a possible discontinuity of your function, or at least its large and oscilating derivations. Is that the case?

[Simon]

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.