I think I solved it....
I got 0.536594(7|8)
My strategy was to use the fact that it was symmetric about x=y.
I took the initial y=... function and rotated it by pi/4 (45 degrees) using the transformation matrix:
Code:
[y'] [cos(pi/4) sin(pi/4)] [y]
[ ] = [ ] [ ]
[x'] [cos(pi/4) -sin(pi/4)] [x]
That allows me to essentially take the integral above the new y' and simply multiply it by 2. I used MATLAB to get it.
Code:
clear
xmin = eval(evalin(symengine, 'numeric::solve(0==y-10*log(y+1)/exp(y+1), y=-0.5..0.5)'));
xmax = eval(evalin(symengine, 'numeric::solve(0==y-10*log(y+1)/exp(y+1), y=0.5..2)'));
error = 1;
stepsize = 3000000;
prevsum = 0;
while (abs(error > 0.000000001))
stepsize = stepsize + 10000;
th = pi/4
x0min = xmin;
x0max = xmax;
y0min = xmin;
y0max = xmax;
x0 = linspace(x0min,x0max,stepsize);
y0 = linspace(y0min,y0max,stepsize);
y = 10*log(x0+1)./exp(x0+1);
x = 10*log(y0+1)./exp(y0+1);
xp = x*cos(th) + y0*sin(th);
yp = y*cos(th) - x0*sin(th);
xp0 = x/cos(th);
plot(x0,y,x,y0,xp,yp);
axis equal
axis square
grid on
sum = 0;
for i=1:stepsize-1
sum = sum + (xp0(i+1)-xp0(i))*yp(i);
end
sum = sum*2;
error = sum - prevsum;
disp('step size');
disp(stepsize);
disp('sum');
disp(sum);
disp('error');
disp(error);
end