Question: An integration issue

I think the worksheet below is enough to define the problem I'm facing.

(question updated by adding a simple 2D case ; the expected value is about 0.67684)

restart

# How can I define CharFunc so that int(CharFunc, x=...) gives me the expected value?
#
# Motivation:
# Let Dom some 2D domain and Env another 2D domain which contains Dom.
# The interior of Dom is defined by a sequence L(x, y) of inequalities.
# The area of Dom can be expressed by the integral over Env of the characteristic function
# CharFunc(x, y) of Dom which returns 1 if (x, y) belongs to Dom and 0 otherwise.
#
# Here is a simple 1D example


CharFunc := proc(x)
  if x::numeric then
    piecewise(is(And(x >= 0, x <= 1)), 1, 0)
  else
    'procname( _passed )'
  end if:
end proc:


Env := -1..2:
plot(CharFunc(x), x=Env, thickness=5);

 

# Expected value = 1

int(CharFunc(x), x=Env);

int(CharFunc(x), x = -1 .. 2)

(1)

# Expected value = 1/2

int(CharFunc(x), x=1/4..3/4);

int(CharFunc(x), x = 1/4 .. 3/4)

(2)


Works with floats

# Expected value = 0.5

int(CharFunc(x), x=0.25..0.75);

.5000000000

(3)

# Expected value = 1

int(CharFunc(x), x=-1.0..2.0);

1.000000000

(4)


2D example

Env := x=0.8..3, y=0..1.3;

x = .8 .. 3, y = 0 .. 1.3

(5)

Dom := proc (x, y) options operator, arrow; And(y <= 1/(1+sinh(2*x)*ln(x)^2), .8 <= x, x <= 3, 0 <= y, y <= 1.1) end proc

proc (x, y) options operator, arrow; And(y <= 1/(1+sinh(2*x)*ln(x)^2), .8 <= x, x <= 3, 0 <= y, y <= 1.1) end proc

(6)

CharFunc := proc(x, y)
  piecewise(is(Dom(x, y)), 1, 0)
end proc:

plot3d(
  [0, 'CharFunc'(x, y)], Env
  , grid=[40, 40]
  , style=surface, color=[gray, blue], transparency=[0, 0]
  , title=typeset('CharFunc'(x, y))
);

 

int(CharFunc(x, y), Env)

0.

(7)
 

 

Download Integration_issue.mw

Thanks in advance

Please Wait...