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.
