How to solve a threedimensional integral of a nested function?

6 次查看(过去 30 天)
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
Mike Hosea 2013-12-16
编辑:Mike Hosea 2013-12-16
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')
  4 个评论
Simon
Simon 2014-1-6
编辑:Simon 2014-1-6
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
Mike Hosea 2014-1-7
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.

请先登录,再进行评论。

更多回答(1 个)

Vaclav Rimal
Vaclav Rimal 2013-12-16
编辑:Vaclav Rimal 2013-12-16
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
  3 个评论
Vaclav Rimal
Vaclav Rimal 2013-12-16
编辑:Vaclav Rimal 2013-12-16
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
Simon 2013-12-18
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.

请先登录,再进行评论。

类别

Help CenterFile Exchange 中查找有关 Numerical Integration and Differentiation 的更多信息

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by