17 views (last 30 days)

Show older comments

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!

Unai San Miguel
on 20 Mar 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).

Torsten
on 20 Mar 2018

Edited: Torsten
on 20 Mar 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.

Torsten
on 20 Mar 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 ?

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

Start Hunting!