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
f=numerically evaluated function using above mesh (so this is a 4D matrix)
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));
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!