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
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).
See Also
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!