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

15 views (last 30 days)

Show older comments

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

##### 5 Comments

Walter Roberson
on 31 Jul 2024

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.

### Accepted Answer

Torsten
on 24 Aug 2024

Moved: Torsten
on 3 Sep 2024

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

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

##### 5 Comments

Torsten
on 26 Aug 2024

Moved: Torsten
on 3 Sep 2024

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.

### More Answers (2)

David Goodmanson
on 2 Aug 2024

Edited: David Goodmanson
on 2 Aug 2024

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.

##### 7 Comments

Walter Roberson
on 16 Aug 2024

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
on 23 Aug 2024

##### 12 Comments

Walter Roberson
on 5 Sep 2024 at 20:34

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
on 5 Sep 2024 at 21:50

Edited: Paul
on 5 Sep 2024 at 21:51

### See Also

### Categories

### Community Treasure Hunt

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

Start Hunting!