double and triple integral for a complex function. I coded but did not get any output

f = (-1) * (x + L/2) ./ ((x + L/2).^2 + (y - y1).^2 + (z - z1).^2).^(3/2);

The goal is to do the double integral of that function over dy1 and dz1 and then the next step is to do the triple integral of the result(double integration result) over dx,dy and dz. I tried these in several ways but it looks like Matlab is unable to to this integration. Any suggestion?

I just want to clarify again. The integration order can not be changed. The integration order should be dy1 dz1 and then dx dy dz. That means first I need to do the double integral (dy1 dz1) of that function and then I need to do the triple integral (dx dy dz) of that double integral result.

I have tried this in MATLAB by using integral2 and integral3 function but nothing works. The 5D integration needs to be solved numerically without changing the integration order.

I have tried this in mathmetica and found out the double integral (dy1 dz1) does not have any closed form solution.

My 3rd step is to try this 5D integral by using integral 2 and then cumetrapze for the triple integral. I can do the first double integral by integral2 function if I define x, y and z as matrix/vectors. I know the limits of x,y and z. So, I can make materices of x, y and z. Then integral2 will retun a results containing scaler numbers. But the problem is now, how I can do the triple integral of this results using cumtrapze. As my double integral answer only contains numbers. How I can convert this number in any form of equation and do the integral with respect to x, y and z. Do I need to write my own integration routine?

Any idea or suggestions?

The limits:

this is for rectangular material layers. here x = length =L = 12e-3, y = width = 3.8 e-3 and z = thickness = tm = 250e-6

y1 = width = 3.8e-3, z1 = thickness of another material layer = tf = 23e-6

then the limit should be y1_max = width/2, y1_min = -width/2, z1_max = tf/2, z1_min = -tf/2,

x_max = length/2, x_min = -length/2, y_max = width/2, y_min = -width/2, z_max = tm/2, z_min = -tm/2

[David Goodmanson]

Hi Orpita, what is the region of integration?

[ORPITA SAHA]

Hi David, this is for rectangular material layers. here x = length = 12e-3, y = width = 3.8 e-3 and z = thickness = tm = 250e-6

y1 = width = 3.8e-3, z1 = thickness of another material layer = tf = 23e-6

then the limit should be y1_max = width/2, y1_min = -width/2, z1_max = tf/2, z1_min = -tf/2,

x_max = length/2, x_min = -length/2, y_max = width/2, y_min = -width/2, z_max = tm/2, z_min = -tm/2

[Paul]

The doc page integral3 shows how to extend functionality to a 4D integral. That can be extended to 6D. That same doc page, in the Tips section, provides a link to the FEX submission integralN that, I believe, can handle up to 6D integrals.

[Walter Roberson]

In theory,

syms x y y1 z z1

syms xl xu yl yu y1l y1u zl zu z1l z1u
f = (-1) * (x + 1/2) ./ ((x + 1/2).^2 + (y - y1).^2 + (z - z1).^2).^(3/2);

step1 =  int(f, y1, y1l, y1u, hold=true)
step2 = int(step1, z1, z1l, z1u, hold=true)
step3 = int(step2, x, xl, xu, hold=true)
step4 = int(step3, y, yl, yu, hold=true)
step5 = int(step4, z, zl, zu, hold=true)

result = release(step5)

In reality, this takes a long time.

[ORPITA SAHA]

Thank you all for the suggestion.

Answer 1

[Torsten]

According to the literature, the Monte Carlo Integration method seems to be the usual approach:

https://iopscience.iop.org/article/10.1088/1757-899X/504/1/012064/pdf

rng("default")
length = 12e-3;
width = 3.8e-3 ;
tm = 250e-6;
y1 = 3.8e-3;
z1 = 23e-6;
y1min = -y1/2;
y1max = y1/2;
z1min = -z1/2;
z1max = z1/2;
xmin = -length/2;
xmax = length/2;
ymin = -width/2;
ymax = width/2;
zmin = -tm/2;
zmax = tm/2;
ny1 = 40;
nz1 = 40;
nx = 40;
ny = 40;
nz = 40;
n = ny1*nz1*nx*ny*nz;
V = y1*z1*length*width*tm;
f = @(y1,z1,x,y,z) (-1) * (x + length/2) ./ ((x + length/2).^2 + (y - y1).^2 + (z - z1).^2).^(3/2);
N = 32;
result = zeros(1,N);
for i = 1:N
  y1 = y1min + rand(ny1,1)*(y1max-y1min);
  z1 = z1min + rand(nz1,1)*(z1max-z1min);
  x = xmin + rand(nx,1)*(xmax-xmin);
  y = ymin + rand(ny,1)*(ymax-ymin);
  z = zmin + rand(nz,1)*(zmax-zmin);
  [Y1,Z1,X,Y,Z] = ndgrid(y1,z1,x,y,z);
  result(i) = sum(f(Y1,Z1,X,Y,Z),'all')*V/n;
end
plot(1:N,cumsum(result)./(1:N))

[ORPITA SAHA]

@Torsten

Thanks for pointing out this way. At N=32, the value of this 5d integral is around -1.6284e-10, which is pretty similar with the value I got by other method (integration order change). But, I think, I can increase the N value to have a more flat (true value) curve with this method.

One more thing, I can not run this code in my installed MATLAB due to out of memory, an error has occured. But, I can successfully ran this code in MATLAB online. Do you recommend to run this code online?

Also, I am not finding the code online when I click the tab 'open in matlab online', . I found earlier. Then I did some minor changes with this code and ran the code. But, now, I can not open this code online. Any suggestion?

[Torsten]

I changed the code from above - it will work now for whatever value of N you choose.

I evaluate f in N packages of size 40^5 grid points and average afterwards.

Of course, computation time scales linearily with N.

[ORPITA SAHA]

Thanks.

[Torsten]

If 40^5 grid points still use too much memory on your personal computer, you can reduce 40 to - say - 20 and increase N accordingly.

Since the value of your integral is that small, integral2 or similar tools are difficult to use since the value of the integral is in the order of the allowed error tolerance. You will have to use very small values for RelTol and AbsTol to get good results, I guess.

[ORPITA SAHA]

Thanks. This method should be appropriate. If you print the codes in the answer box here, I can accept it.

Answer 2

[David Goodmanson]

Hi Orpita,

Here is a way to at least reduce the dimension of the integral from 5d to 3d. First of all, the quantity (x+1/2) seems like an unusual choice since the interval of integration for x is only

(-1/2)*12e-3 <= x <= (1/2)*12e-3 

It's good if it is 1/2 though since that way the integrand can never have any near-singularities due to a tiny denominator. Here I will just use x0, x0 = 1/2 in this case, and also use xa and xb as the min and max limits for the x integration and similarly with ya,yb etc. for the others.

To cut to the chase, the 3d integral is

Int g(y1,z,z1) dy1 dz dz1                                                                (2)
where
g(y1,z,z1) = asinh((yb-y1)/Ab) -asinh((ya-y1)/Ab) -asinh((yb-y1)/Aa) + asinh((ya-y1)/Aa)
and 
Ab = ( (xb+x0).^2 + (z-z1).^2 )^(1/2) 
Aa = ( (xa+x0).^2 + (z-z1).^2 )^(1/2)

Boring details: The initial integral is

Int  (-1) * (x+x0)./( (x+x0).^2 + (y-y1).^2 + (z-z1).^2 ).^(3/2)  dx dy dy1 dz dz1

The integrand is a perfect differential in x so we can immediately do the x integration to obtain the 4d integral

Int    [ 1/ ((xb+x0).^2 + (y-y1).^2 + (z-z1).^2).^(1/2) ...
        -1/ ((xa+x0).^2 + (y-y1).^2 + (z-z1).^2).^(1/2) ]        dy dy1 dz dz1     (1)

Two basically identical integrands here. Consider the top one, and one can choose which variable to integrate next. Lets say y. Denote

Ab^2 = (xb+x0).^2 + (z-z1).^2  

and make the substitution

(y-y1) = Ab sinh(u)   d(y-y1) = Ab cosh(u)  du      % here y1 is treated as a constant.

Plug that in, obtain

Int Ab cosh(u) du  / (Ab^2 + Ab^2 sinh^2(u))^(1/2) = Int Abcosh(u)/Abcosh(u)du = u
% and still dy1 dz dz1

now u = asinh( (y-y1) /Ab) and this is evaluated at the limits of y to obtain

asinh((yb-y1)/Ab) - asinh((ya-y1)/Ab)

This goes the same for the lower integral in (1) so the final result is as shown in (2).

As a side note, it is not a good idea to have a variable called 'length' since that is the name of a Matlab function.

[ORPITA SAHA]

Hi David, Thank you.I think there is a misunderstanding about the main equation. The dinominator of the equation is (x + L/2), its not (x+1/2). L = length = 12e-3. I used the lower case of L, that's why this misunderstanding happened. Same thing applies to the numenator. The first term of the numinator is ((x +L/2)^2 +........

Now, how would you solve?

I thought a way. Like first using integral2 function for dy1 and dz1. As itegral2 results need to be scaler, I have to use marix of x, y and z. I mean I know the limits of them. So I can make matrices of x, y and z. And then I will use "cumtrapze" for integrating that with x, y and z. Just got this idea did not implement that yet. Do you have any thoughts on that?

Previously, I have used wolfram mathematica to solve that integral and found out that integral does not have any closed form solution. Also, as my L is so small, the numerical integration of mathematica does not provide an appropriate results. So, this numerical integration needs to be done by MATLAB.

[Torsten]

x0 can also be set to L/2 instead of 1/2. So as far as I can see, @David Goodmanson 's analysis is still valid.

@David Goodmanson 's analysis leaves you with a 3d-integral over dy1 dz dz1 that you could try to evaluate using "integral3".

[ORPITA SAHA]

Thanks @David Goodmanson and @Torsten.

[David Goodmanson]

Hi Orpita,

I should have spotted the small l, but for unclear-font reasons like this, I myself never use a variable called small l, and usually (unless it's clear from the context) any variable name that includes a small l. Notice how in this very message with its sans serif type, you can't tell the difference at all berween small l and capital I.

It appears that you have

X = (x+L/2)   with  -L/2 <= x <= L/2

so it's possible for X to take on the value 0, is that correct? If so, then when y=y1 and z=z1 there is a singularity in

f = (-1) * (x+L/2) ./ ((x+L/2).^2 + (y-y1).^2 + (z-z1).^2).^(3/2);

due to a zero in the denominator. But it is not a strong singularity and the function otherwise is smooth and nonoscillatory so I like the idea of using cumtrapz in three dimensions, especially if you change the boundary conditions very slightly so that you can make the grid spacings for e.g. z and z1 to be noncommensurate, so that (z-z1) can't be zero on the grid.

If neither of your analytic integrals is done in x I think you are going to have a difficult time finding an analytic exppression for the integration of any two of y,y1,z,z1.

[ORPITA SAHA]

Hi David,

I did the 5d integration and got a value with the method you have shown earlier (first dx and dy, this gives a closed form solution, and then I used this equation and solved it by integral3 function). But, I have a question here. In this method, I am besically changing the order of my integration. I mean I am doing the integration first dx, then dy and then dy1dzdz1. The actual integration order should be dy1dz1 and then dxdydz. That's why, I suspect the value I got with the method (frrst dx, 2nd dy and then dy1dzdz1) is correct or not? So far I know, if the function has singularity, I can not change the order of integration. And, my function has singularity.

@David Goodmanson what's your thoughts on this?

[Walter Roberson]

By one of the fundamental theorems of integration, the order of integration does not matter, except for cases where the limits depend upon one of the variables being integrated over.

If the function has a singularity, you cannot integrate it, unless you are working with something like the Cauchy Principle Value

[ORPITA SAHA]

@Walter Roberson Thanks. But, I think, If I change the order of integration this fuction will not behave appropriately. Because, without changing the order, the integral does not have any closed form solution. If I change the order of integration it has a closed form solution. It is suspicious. Also, I am not getting my arropriate value at the end in this method. So, I can not change the order of integration. This function contains singularity.

Answer 3

[ORPITA SAHA]

I just want to clarify again. The integration order can not be changed. The integration order should be dy1 dz1 and then dx dy dz. That means first I need to do the double integral (dy1 dz1) of that function and then I need to do the triple integral (dx dy dz) of that double integral result.

I have tried this in MATLAB by using integral2 and integral3 function but nothing works. The 5D integration needs to be solved numerically without changing the integration order.

I have tried this in mathmetica and found out the double integral (dy1 dz1) does not have any closed form solution.

My 3rd step is to try this 5D integral by using integral 2 and then cumetrapze for the triple integral. I can do the first double integral by integral2 function if I define x, y and z as matrix/vectors. I know the limits of x,y and z. So, I can make materices of x, y and z. Then integral2 will retun a results containing scaler numbers. But the problem is now, how I can do the triple integral of this results using cumtrapze. As my double integral answer only contains numbers. How I can convert this number in any form of equation and do the integral with respect to x, y and z. Do I need to write my own integration routine?

Any idea or suggestions?

The limits:

this is for rectangular material layers. here x = length = 12e-3, y = width = 3.8 e-3 and z = thickness = tm = 250e-6

y1 = width = 3.8e-3, z1 = thickness of another material layer = tf = 23e-6

then the limit should be y1_max = width/2, y1_min = -width/2, z1_max = tf/2, z1_min = -tf/2,

x_max = length/2, x_min = -length/2, y_max = width/2, y_min = -width/2, z_max = tm/2, z_min = -tm/2

[Walter Roberson]

Your integrals all have constant bounds. According to basic theorems https://en.wikipedia.org/wiki/Order_of_integration_(calculus)#Basic_theorems the order of integration is irrelevant, provided that the functions are continuous and integrable.

Because, without changing the order, the integral does not have any closed form solution.

The int() function essentially performs pattern matching to try to determine analytic integrals. Sometimes, integrating by one variable first does not happen to match any of the stored patterns, but integrating by a different variable does happen to match one of the stored patterns. Therefore, int() over one variable might fail to find an analytic solution whereas int() over a different variable might find an analytic solution. This does not mean that the function has no analytic integral if done in one order instead of the other: it merely means that int() was not powerful enough to find the analyatic integral.

[ORPITA SAHA]

yes, That's true. But, my function conatins sigularity when x = L/2, y =y1 and z = z1. If a function has singularities, it will not give you the same results if you change the order of integration.

[Torsten]

@ORPITA SAHA

I didn't spend time checking if your 5d integral exists. But if you want to try what you outlined above, I suggest the following loop. But you will have to check whether the final result stabilizes for different values of nx, ny and nz.

xmin = -length/2;
xmax = length/2;
ymin = -width/2;
ymax = width/2;
zmin = -tm/2;
zmax = tm/2;
nx = 10;   % Number of intervals in x-direction
ny = 10;   % Number of intervals in y-direction
nz = 10;   % Number of intervals in z-direction
dx = (xmax-xmin)/nx;   % Spacing in x-direction
dy = (ymax-ymin)/ny;   % Spacing in y-direction
dz = (zmax-zmin)/nz;   % Spacing in z-direction
for ix = 1:nx
    xmid = dx/2 + (ix-1)*dx;
    for iy = 1:ny
        ymid = dy/2 + (iy-1)*dy;
        for iz = 1:nz
            zmid = dz/2 + (iz-1)*dz;
            % Compute 2d-integral over y1 and z1 in (xmid,ymid,zmid)
            two_dimensional_integral(ix,iy,iz) = ...;
        end
    end
end
% Approximate 5d-integral 
fifth_dimensional_integral = sum(two_dimensional_integral(:))*dx*dy*dz

[Walter Roberson]

integral() over a singularity might return anything.

If you are after the Cauchy Principle Value, then you need to split the interval and add the results of separate integrations with the singularities at the endpoints.

[Torsten]

Before starting any numerical computations, it must be established that the 5d - integral exists (despite the singularities). If the integral has a physical meaning, asking for advice in a physics forum or at least in a mathematics forum makes sense in my opinion.

[ORPITA SAHA]

@Torsten @Walter Roberson

Thanks. The integral is taken from a published Journal. It is one of the model equations for a Magnetoelectric laminate.

[Torsten]

Is the integral evaluated in the journal ? Can you provide a link ?

[ORPITA SAHA]

No the integral is not evaluated there.

https://www.sciencedirect.com/science/article/pii/S0263822314004991

Please see the above link. See the equation number (2), (3) and (4). If you put the equations of Hf3(+) and Hf3(-) to the equation (2), it forms 5D integral.

I am also attaching the journal here.

[Paul]

@Torsten

From this comment: "Before starting any numerical computations, it must be established that the 5d - integral exists."

Can you expand on how one might establish the existence of the integral, or any integral in general, and if Matlab can be used in that process?

[Torsten]

Can you expand on how one might establish the existence of the integral, or any integral in general, and if Matlab can be used in that process?

I don't think that MATLAB can be used for this.

Usually, theoretical results from Lebesgue integration theory are used to establish that functions are integrable, e.g. the Dominated Convergence Theorem:

https://en.wikipedia.org/wiki/Dominated_convergence_theorem

But in the present case, I can't come up with a way to proof the existence of the 5-fold integral. According to the publications in this field, it seems that it's obvious :-)

[Walter Roberson]

integral() over a singularity might return anything

If you are using the integral*() family of functions, you must split the calculations at the boundary conditions.

[Paul]

@Torsten

My recollection of the basics is that, for the 1D case, a sufficient condition for the integral to exist is the limits of integration (a and b) are finite and the integrand f(x) is continuous on the closed interval [a,b]. Now that I think of it, I'm not even sure if Matlab can be used to verifty that a function to be integrated, either with int or integral, is continuous. There might be easy cases on the symbolic side to show that f(x) is discontinuous (e.g., if f(x) includes a heaviside).