How to integrate when the integrand is not a built-in function?

Posted on 2008-06-22 14:58 By dpaddy (

15)

Categories:

How do I ...? with Maple (Student)

I have a nice function whose graph is approximately quadratic; the plot shows the area under it should be around 1.6

However, it seems Maple thinks the area is 0.1465

Presumably there is something simple I am doing wrong, but I can't see it.

What follows may seem strange -- it is certainly strange to me -- but I need to be able to make integration work with my functions as apposed to the built-in functions one sees in examples...

with(Student[Calculus1]);

Digits:= 64;

erfy := (x) -> piecewise(x < -5, -1+(-1/Pi^(1/2)/x+1/2/Pi^(1/2)/x^3-3/4/Pi^(1/2)/x^5)/exp(x^2), And(-5 <= x, x <= 5), erf(x), x > 5, (1+(-1/Pi^(1/2)/x+1/2/Pi^(1/2)/x^3-3/4/Pi^(1/2)/x^5)/exp(x^2)));

h:= (u,a,b,g,d) -> exp(-(1/d)^2)*d*g*exp(-(u-a)^2/b^2)*exp(((d-u)/(d*u))^2)*(1-erfy((d-u)/(d*u)))/u;

core := proc(n,t0,t1,a,b,g,d)

local x,y,f;

x := 2*t0/(d*g);

y := 2*t1/(d*g);

f := unapply(h(u,a,b,g,d), u);

[ApproximateInt(f, x..y, method = simpson),x,y];

end proc;

core(3,.5, 2.1, .5, 3, 3.14, 2.4);

MAPLE responds with an integral of about 0.1465...

HOWEVER:

show := proc(n,t0,t1,a,b,g,d)

local x,y,f;

x := 2*t0/(d*g);

y := 2*t1/(d*g);

f := unapply(h(u,a,b,g,d), u);

plot(f, x..y);

end proc;

show(3,.5, 2.1, .5, 3, 3.14, 2.4);

MAPLE responds with the nice quadratic graph whose area should be about 1.6...

Is it possible that Maple can be so wrong?!? More likely would be that I have no clue as to how one somehow manages to cajole Maple into preforming an integral.

Please help.

1.6366587497074052272870953240180323794215383 ...

Comment on 2008-06-22 16:10 from Axel Vogt ( Maple Rating 4

994)

Have not stepped through the problem in detail, but you do not need
to approximate the error function on your own (it makes the function
non-analytic) and there is no need to use Simpson (playing with your
approach I see a 'discontinuity' at ~ 1.6 for your task, but do not
want to look closer).

After converting to rationals you want the following:

Int(169858969913624935089195720356683/26813033068808184917956157510272*
(1-erf(5/12*(12/5-u)/u))/u*exp(-1/9*(u-1/2)^2+25/144*(12/5-u)^2/u^2),
u = 125/942 .. 175/314);

Using simplify(%); expand(%); combine(%,exp); normal(%); evalf(%);
gives the value in the title (astonishingly not 64 digits ...)

ApproximateInt

Comment on 2008-06-22 16:25 from Robert Israel ( Maple Rating 4

1555)

ApproximateInt works correctly if you use an expression instead of a function, i.e.

ApproximateInt(f(t),t=x..y,method=simpson)

However, the syntax ApproximateInt(f, x..y) does not work. In fact, the result does not depend on f at all.

> ApproximateInt(f, a .. b);

1/2*b^2-1/2*a^2

[ The <maple> tag is acting up again: that should be 1/2*b^2 - 1/2*a^2 ]

What seems to be happening is that Maple interprets this input as

ApproximateInt(f, f=a..b)

To be fair, the help page ?ApproximateInt does not mention this syntax at all, so ApproximateInt(f, x..y) is an error on the part of the user. However, this is definitely something that should be caught (either returning the same result as ApproximateInt(f(t), t=x..y) or an error message). I'm submitting an SCR.

evalf/Int

Comment on 2008-06-22 16:42 from alec ( Maple Rating 4

1608)

The more or less standard way is to use evalf/Int. Somehing like

core1 := (t0,t1,a,b,g,d,eps)->
evalf(Int(h(u,a,b,g,d), u=2*t0/(d*g)..2*t1/(d*g),epsilon=eps)):

core1(.5, 2.1, .5, 3, 3.14, 2.4, 10^(-32));
 1.636668036549567312072601247206393398244078914859354065223268055

It is interesting that the answer is different than in Axel Vogt's post above.

Another bug related to this is

Digits:=10:
core2 := (t0,t1,a,b,g,d,dig)->
evalf(Int(h(u,a,b,g,d), u=2*t0/(d*g)..2*t1/(d*g),digits=dig)):

core2(.5, 2.1, .5, 3, 3.14, 2.4, 6);
                                       14
                            0.304448 10

(But I am not submitting the SCR.)

Alec

looking closer

Comment on 2008-06-23 14:38 from Axel Vogt ( Maple Rating 4

994)

May be, because I use erf instead of his approximation 'erfy'
and was working with rationals ...

I think the numerical problem can be understood: He uses a term
1 - erf(  (((d-u)/(d*u))) ) in a range 0.13 ... 0.55 where the
inner term is (2.4-u)/u/2.4 and plotting let it reach infinity
towards u ~ 0, it is falling from ~ 7 to 0 in that range.

My ~1.6 above is a typo - it is ~ 0.16. This for the one point
is, where his 'erfy' switches. The other is: then one gets some
precision & cancellation  problems by 1 - erf(_large_) ~ 0 and
it is better to use the complementary erfc.

Additionally that is multiplied by exp( innerTerm^2 ) achieving
exp(7^2) ~ 10^22 and multiplied a Gaussian exp(-(u-a)^2/b^2),
which is ~ 1 in the range and erfc falls down to 10^(-24).

Now Maple is going to evaluate that functions ... and the exp
almost explodes at the left boundary with flat erfc.

After that feeling says, that for integration it must be handled
as a numerical singularity.

Note that a plot will not care to much for it (using precision
high enough), but an integrator may got trapped.


For the task exp(-(1/d)^2)*d*g is a constant, which should be
done outside the integral.

Note that the parameter n=3 in the function 'core' is *not used*.

Mills' ratio

Comment on 2008-06-23 15:13 from Axel Vogt ( Maple Rating 4

994)

After some 3rd thoughts ... the term exp(((d-u)/(d*u))^2)*(1-erf((d-u)/(d*u)))
is Mills' ratio I think, for which *good* numerics are available (I did that
in the code of my blog posts on exact computations for Black-Scholes and the
computation of implied volatility).

So the integrand is like Gaussian * Mills ratio and for proper integration now
I would code the ratio as separate function (it is a friendly one) and the
task should be done in usual double precision, i.e. with 14 Digits.

Thanks to all!

Comment on 2008-06-22 17:48 from dpaddy (

15)

Thanks to all!