multidimensional numerical integration without large matrices

9 visualizaciones (últimos 30 días)
kristoffer svendsen
kristoffer svendsen el 19 de Mzo. de 2018
Comentada: Torsten el 20 de Mzo. de 2018
I have a function of 4 variables which neeeds to be integrated over a 3D space. The integral does not have an analytical solution so it needs to be done numerically. The function is also very sensitive and easily becomes unstable if the step size is too small. The way I've done it previously is as
[X1,X2,X3,X4]=meshgrid(x1,x2,x3,x4)
f=numerically evaluated function using above mesh (so this is a 4D matrix)
trapz(X1,f,3)
trapz(X1,f,2)
trapz(X1,f,1)
This however, bottlenecks very fast when I need to increase the evaluation points to make it stable (runs out of memory). So my question is if there is another way to do this, without having 4D matrices.
I tried something like
g1=@(x1,x2,x3) integral(@(x4) f(x1,x2,x3,x4,x4(1),x4(end));
g2=@(x1,x2) integral(@(x3) g1(x1,x2,x3,x3(1),x3(end));
g3=@(x1) integral(@(x1) g2(x1,x2,x3,x4, x2(1),x2(end));
g3(x1)
which does not work (gets matrix dimension mismatch error despite correctly vectorized function), as well as
g=@(x1) integral3(@(x2,x3,x4) F(x1,x2,x3,x4),x2(1),x2(end),x3(1),x3(end),x4(1),x4(end))
which fails the global error test, returning NaN.
The function in question is
F=@(x1,x2,x3,x4) sqrt(c./(1i.*x4.*a.*b)) .* exp((-1i.*pi./(x4.*b)).* (2.*x3.*x1+b./a.*x2-c./a.*x3.^2-a./c.*(x1+b./a.*x2).^2) ).*heaviside(abs(x3)-d) ;
a,b,c,d are constants
Would highly appriciate any suggestions!
  2 comentarios
Walter Roberson
Walter Roberson el 19 de Mzo. de 2018
"which fails the global error test"
? Could you explain that a bit more?
kristoffer svendsen
kristoffer svendsen el 20 de Mzo. de 2018
It states "Warning: Reached the maximum number of function evaluations (10000). The result fails the global error test."
By decreasing the tolerances this is solved but the integration is incorrect as it shows only noise.

Iniciar sesión para comentar.

Respuestas (3)

Walter Roberson
Walter Roberson el 19 de Mzo. de 2018
  1 comentario
kristoffer svendsen
kristoffer svendsen el 20 de Mzo. de 2018
Tried this but this integrates over all variables, I want to integrate over 3 out of 4 variables to be left with either a vector or another function F @(x1). Maybe I misunderstood this script but I could not get it to work this way.

Iniciar sesión para comentar.


Unai San Miguel
Unai San Miguel el 20 de Mzo. de 2018
I don't like meshgrid even for 2 variables, so I would use ndgrid instead
[X1, X2, X3, X4] = ndgrid(x1, x2, x3, x4);
This way X1(i, j, k, l) = x1(i), X2(i, j, k, l) = x2(j) ans so on. Your array f would be of size [n1, n2, n3, n4], being n1 = length(x1), .... Then you could use trapz but with the lowercase-single dimensional variables.
g(x1, x2, x3) = int(f, dx4) ~ |trapz(x4, f, 4)|
h(x1, x2) = int(g, dx3) ~ |trapz(x3, trapz(x4, f, 4), 3)|, ...
So your final function would be
g_1 = trapz(x2, trapz(x3, trapz(x4, f, 4), 3), 2);
an array of size n1. I haven't done integrals of more than 2 variables, but if you can handle the 4D f array this looks doable (you are always reducing the size of the arrays).
  1 comentario
kristoffer svendsen
kristoffer svendsen el 20 de Mzo. de 2018
Editada: kristoffer svendsen el 20 de Mzo. de 2018
My apologies, yes I do indeed use ndgrid and not meshgrid, as meshgrid does not accept more than 3 dimensions. The 4D array however takes up about 60 Gb of memory and can not be stored for computation.

Iniciar sesión para comentar.


Torsten
Torsten el 20 de Mzo. de 2018
Editada: Torsten el 20 de Mzo. de 2018
x2min = ...;
x2max = ...;
x3min = ...;
x3max = ...;
x4min = ...;
x4max = ...;
a = ...;
b = ...;
c = ...;
d = ...;
F=@(x1,x2,x3,x4) sqrt(c./(1i.*x4.*a.*b)) .* exp((-1i.*pi./(x4.*b)).* (2.*x3.*x1+b./a.*x2-c./a.*x3.^2-a./c.*(x1+b./a.*x2).^2) ).*heaviside(abs(x3)-d) ;
g=@(x1) integral3(@(x2,x3,x4) F(x1,x2,x3,x4),x2min,x2max,x3min,x3max,x4min,x4max);
g(2.5)
does not work ?
Which values do you use for the "..." indicated constants ?
Best wishes
Torsten.
  2 comentarios
kristoffer svendsen
kristoffer svendsen el 20 de Mzo. de 2018
Editada: kristoffer svendsen el 20 de Mzo. de 2018
This works if the 'AbsTol' and 'RelTol' are decreased (to roughly 1e-3) but at this point the integration does not converge to the correct solution and the result is not as it should be.
The constants varies (hehe) but an example would be
c=0.15; b=0.1; a=0.05; d=25e-6;
x2min=-7.5e-5;
x2max=7.5e-5;
x3min=x2min;
x3max=x2max;
x4min=1.2e-10;
x4max=1.2e-9;
Torsten
Torsten el 20 de Mzo. de 2018
What if you split the integral into three integrals
x3min <= x3 <= -d
-d <= x3 <= d
d <= x3 <= x3max
and remove the heaviside term in F ?

Iniciar sesión para comentar.

Categorías

Más información sobre Numerical Integration and Differentiation en Help Center y File Exchange.

Community Treasure Hunt

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

Start Hunting!

Translated by