4D integral with exponential value

Hi,

I am trying to calculate this integral using indefinite integral functions.

 1   1    /   1   1
    /   /    |   /   /
   |   |     |  |   |                                     1
   |   |  exp|  |   |  ----------------------------------------------------------------------- dx dy
  /   /      | /   /   /                1                     \ /        2   /         1 \2 \2
    0   0    |   0   0 | ------------------------------ + 160 | | (w - z)  + | y - x + - |  |
             |         | /        2   /         1 \2 \2       | \            \         2 /  /
             |         | | (w - z)  + | y - x + - |  |        |
             \         \ \            \         2 /  /        /

     \
     |
     |
     | dz dw
     |
     |
     |
     |
     /

It is 4D integral. When I try to use indefinite integral function int Matlab cant calculate it.

int(int(exp(int(int(1/((1/((w - z)^2 + (y - x + 1/2)^2)^2 + 160)*((w - z)^2 + (y - x + 1/2)^2)^2), x, 0, 1), y, 0, 1)), z, 0, 1), w, 0, 1)

When I try to use integral2 function to calculate numerically, the z and w parameters are obstacles.

I cannot use 4D user defined integral function because there exp{} block between integrals.

Do you know how to solve that integral?

Answer 1

[Torsten]

Brute-force method:

delta=0.01;
w=0:delta:1;
z=0:delta:1;
x=0:delta:1;
y=0:delta:1;
func=@(x,y,w,z)1/(1+160*(((w-z)^2+(y-x+0.5)^2)^2));
for i=1:numel(w)
  wloc = w(i);
  for j=1:numel(z)
    zloc = z(j);
    f = 0;
    for k=1:numel(x)-1
      xloc_l = x(k);
      xloc_r = x(k+1);
      for l=1:numel(y)-1
        yloc_l = y(l);
        yloc_r = y(l+1);
        f = f + 0.25*(func(xloc_l,yloc_l,wloc,zloc)+func(xloc_l,yloc_r,wloc,zloc)+func(xloc_r,yloc_l,wloc,zloc)+func(xloc_r,yloc_r,wloc,zloc))*delta^2;
      end
    end
    fxy(i,j) = exp(f);
  end 
end 
integral = 0.0;
for i=1:numel(w)-1
  for j=1:numel(z)-1
    integral = integral + 0.25*(fxy(i,j)+fxy(i,j+1)+fxy(i+1,j)+fxy(i+1,j+1))*delta^2;
  end 
end

"integral" should be the value of the integral you are looking for.

You may want to choose a smaller value for "delta" to get a better approximation.

With delta=1/300 I get a value of 1,166618 for your integral.

Best wishes

Torsten.

Answer 2

[Mümin ÖZPOLAT]

Thanks !

Seems like it is an implementation of Riemann Sum.

You use delta^2 because the integration is 2 layered and multiplicate it by 0.25 because you get average of 4 different result.

[Torsten]

It's the two-dimensional trapezoidal rule, applied twice: once to the double inner integral, then to the double outer intergral.

Best wishes

Torsten.