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

Question

ORPITA SAHA 2024-7-30

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2141621-double-and-triple-integral-for-a-complex-function-i-coded-but-did-not-get-any-output

编辑： Paul 2024-9-5，21:51

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 个评论
显示 3更早的评论隐藏 3更早的评论

Walter Roberson 2024-7-31

在 MATLAB Online 中打开

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 2024-8-1

Thank you all for the suggestion.

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

Torsten 2024-8-24

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2141621-double-and-triple-integral-for-a-complex-function-i-coded-but-did-not-get-any-output#answer_1510124

移动：Torsten 2024-9-3，20:45

在 MATLAB Online 中打开

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

5 个评论
显示 3更早的评论隐藏 3更早的评论

Torsten 2024-8-26

移动：Torsten 2024-9-3，20:49

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 2024-9-3，19:05

移动：Torsten 2024-9-3，20:49

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

请先登录，再进行评论。

Answer 2

David Goodmanson 2024-8-2

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2141621-double-and-triple-integral-for-a-complex-function-i-coded-but-did-not-get-any-output#answer_1493621

编辑：David Goodmanson 2024-8-2

在 MATLAB Online 中打开

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 个评论
显示 5更早的评论隐藏 5更早的评论

ORPITA SAHA 2024-8-2

编辑：ORPITA SAHA 2024-8-2

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.

David Goodmanson 2024-8-3

编辑：David Goodmanson 2024-8-23

在 MATLAB Online 中打开

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.

Walter Roberson 2024-8-16

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 2024-8-22

@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 2024-8-23

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2141621-double-and-triple-integral-for-a-complex-function-i-coded-but-did-not-get-any-output#answer_1504399

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

12 个评论
显示 10更早的评论隐藏 10更早的评论

Walter Roberson 2024-9-5，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 2024-9-5，21:50

编辑：Paul 2024-9-5，21:51

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

请先登录，再进行评论。

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

5 个评论
显示 3更早的评论隐藏 3更早的评论

采纳的回答

5 个评论
显示 3更早的评论隐藏 3更早的评论

更多回答（2 个）

7 个评论
显示 5更早的评论隐藏 5更早的评论

12 个评论
显示 10更早的评论隐藏 10更早的评论

另请参阅

类别

标签

产品

版本

Community Treasure Hunt

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

5 个评论 显示 3更早的评论隐藏 3更早的评论

采纳的回答

5 个评论 显示 3更早的评论隐藏 3更早的评论

更多回答（2 个）

7 个评论 显示 5更早的评论隐藏 5更早的评论

12 个评论 显示 10更早的评论隐藏 10更早的评论

另请参阅

类别

标签

产品

版本

Community Treasure Hunt

WeChat

5 个评论
显示 3更早的评论隐藏 3更早的评论

5 个评论
显示 3更早的评论隐藏 3更早的评论

7 个评论
显示 5更早的评论隐藏 5更早的评论

12 个评论
显示 10更早的评论隐藏 10更早的评论