On April 14, 2012, Barton Willis wrote:
-----------------------
>bug?
>
> (%i26) mydefint2(1,[y,-x,x],[x,-1,1]);
> (%o26) res79
>
I have provided for such a user error in the new version.
>Also
>
> separate(l,x) := map('listify, partition_set(setify(l), lambda([s],
> freeof(x,s))))
>
> might useful--I think your use of length causes trouble:
>
> (%i28) separate([x,y],x);
> length: argument cannot be a symbol; found x
---------------------------------
The purpose of the function separate was not well described.
It's use is to filter the existing list facts(), each element of which
will (probably?) have the form x < a, or x > 0, or a < b, for
example. I want to ignore elements like a<b and x<a, but
forget elements like x < 0. I have yet to write a filter function
which will ignore an element like x < a.
The present working version does not call separate, and hence
is not assured to be reliable unless none of the global list
facts() involves the outer integral variable (here called x).
The calling syntax mydefint2(expr, [x,x1,x2], [y,y1,y2]) is the
same as adopted by mathematica, namely the y integral is
the inner integral, whose limits can depend on x, and the
x integration is the outer integration, whose limits should
not depend on y.
--------------------------------------
On April 15, Barton Willis wrote:
--------------------------------------
>The integrand (first argument to mydefin1) is automatically simplified in
>the current context; thus
>
> (%i2) (assume(x < 0), mydefint1(abs(x),x,0,1));
> (%o2) -1/2
>
>No amount of quoting changes this. If this bothers you
> (first simplify in global context and second in the
>context defined by the definite integral), you could specify
>the integrand as a lambda form, or you could write
>your code in Common Lisp as a defmspec function.
---------------------------------------
I agree that I need to write mydefint1 and mydefint2 as
defmspec lisp functions.
The present version (which assumes facts() has no
elements involving x, as described above in more
detail) has been modularized to cover two different
cases. If the outer integral limits involve unbound
parameters, the code makes assumptions and informs
the user.
This looks clunky, but this code is not being
developed as a standalone utility, but rather is intended
to only be called by the "nintegrate" utility I have
been calling nint, which is designed to produce a
floating point number, and hence the outer integral
limits will always be bound to numbers.
There are numerical integrals which quadpack cannot
do, but integrate can do, so unless there are other problems
such as multivalued expressions, etc., it is useful to first
try integrate and see if it succeeds.
In so doing, I want to minimize the number of variable
questions which can easily be answered by looking at
the limits requested.
The present version of mydefint2 then produces:
-----------------------------------
(%i1) load(mydefint);
(%o1) "c:/work2/mydefint.mac"
(%i2) mydefint2(1,[x,0,1],[y,0,1]);
(%o2) 1
(%i3) mydefint2(1,[x,0,a],[y,0,b]);
Since the outer integral limits involve unbound parameter(s),
the code assumes that
0 <= x <= a
(%o3) a*b
(%i4) mydefint2(x*y,[x,0,1],[y,0,1]);
(%o4) 1/4
(%i5) mydefint2 (log(x+y),[x,0,1],[y,0,1]);
(%o5) (4*log(2)-3)/2
(%i6) mydefint2 (1/sqrt(x+y),[x,0,1],[y,0,1]);
(%o6) (2^(7/2)-8)/3
(%i7) mydefint2(sqrt(x^2 - 2*x*y + y^2),[x,-1,1],[y,-1,1]);
(%o7) 8/3
(%i8) mydefint2(sqrt(x^2 - 2*x*y + y^2),[x,-1,1],[y,-x,x]);
(%o8) 0
(%i9) mydefint2(1,[x,-y,y],[y,-1,1]);
the outer integral limits should not depend
on your chosen inner integral variable of integration
(%o9) false
-------------------------------------
and the present code is:
----------------------------------------
/* mydefint.mac */
/* future work
load("c:/work2/mydefint.lisp")$
*/
/* separate(aL,x) returns a list of the elements
of aL which involve x .
to do: need to design a filter which looks for
symbolic limits like x < a and ignores
them
separate(a11L,xx11) :=
block ( [t11,t11L:[]],
for t11 in a11L do
if (part (t11,1) = xx11 or part (t11,2) = xx11) then
t11L : cons (t11,t11L),
t11L)$
*/
/* mydefint1 is only reliable if facts() = []
or if facts() has no entries involving
the variable of integration .
to do: write as a defmspec lisp function.
*/
mydefint1 (e78, x78, a78, b78) :=
block( [domain : complex, res1, loc_assume,e781],
loc_assume : apply ('assume, [x78 > min(a78,b78), x78 <
max(a78,b78)] ),
e781 : expand (e78,0,0),
res1 : ratsimp (integrate (e781, x78, a78, b78 ) ),
apply ('forget, loc_assume),
res1)$
/* mydefint2_symb_lim makes assumptions about
ordering of outer integral variable x relative
to x2 and x1 if qxL = [x,x1,x2], say, and
informs user.
This code is only reliable if global facts()
= [] or entries don't involve qxL[1].
Calls mydefint1.
to do: write as a defmspec lisp function. */
mydefint2_symb_lim(qe,qxL,qyL) :=
block ( [locx_assume,qe79,res79],
/* this is the case in which the outer integral limits
involve unbound parameters */
print (" Since the outer integral limits involve unbound parameter(s),
the code assumes that ",qxL[2]," <= ",qxL[1]," <= ",qxL[3]),
locx_assume : apply ('assume, [qxL[1] > qxL[2],
qxL[1] < qxL[3]] ),
qe79 : expand (qe,0,0),
res79 : apply ('integrate,
flatten ( [apply ('mydefint1,flatten ( [qe79,qyL])), qxL])),
apply ('forget, locx_assume),
ratsimp (res79))$
/* mydefint2_numer_lim is only reliable if global facts()
= [] or entries don't involve qxL[1].
Calls mydefint1.
to do: write as a defmspec lisp function. */
mydefint2_numer_lim (qe,qxL,qyL) :=
block ( [locx_assume, qe79,
qxmin,qxmax,res79,res791,res792,xsign],
/* this is the case in which the outer integral limits are
numbers, or involve %pi or %e etc which will reduce to
numbers when they appear as args of min and max and
compared to each other or to numbers :
(%i18) max(1,%pi);
(%o18) %pi */
qxmin : min (qxL[2], qxL[3]),
qxmax : max (qxL[2], qxL[3]),
if qxmin = qxL[3] then xsign : -1 else xsign : 1,
/* case x doesn't change sign */
if (qxmin >= 0 or qxmax <= 0) then
(locx_assume : apply ('assume, [qxL[1] > qxmin,
qxL[1] < qxmax] ),
qe79 : expand (qe,0,0),
res79 : apply ('integrate,
[apply ('mydefint1,flatten ( [qe79,qyL])), qxL[1], qxmin,
qxmax]),
res79 : xsign*ratsimp (res79),
apply ('forget, locx_assume))
else /* case x does change sign: split into
two integrals */
( locx_assume : apply ('assume, [qxL[1] > qxmin,qxL[1] < 0] ),
qe79 : expand (qe,0,0),
res791 : apply ('integrate,
[apply ('mydefint1, flatten ([qe79,qyL])),qxL[1],qxmin,0]),
apply ('forget, locx_assume),
locx_assume : apply ('assume, [qxL[1] > 0,qxL[1] < qxmax] ),
qe79 : expand (qe,0,0),
res792 : apply ('integrate,
[apply ('mydefint1,flatten ([qe79,qyL])),qxL[1],0,qxmax]),
res79 : xsign*ratsimp (res791 + res792),
apply ('forget,locx_assume)),
res79)$
/** mydefint2 is only reliable if global facts()
= [] or entries don't involve qxL[1].
Calls either mydefint2_numer_lim
or mydefint2_symb_lim.
to do: write as a defmspec lisp function. */
mydefint2 (pqe,pqxL,pqyL) :=
block ( if member (pqyL[1], listofvars(pqxL)) then
(print ("the outer integral limits should not depend
on your chosen inner integral variable of integration"),
return (false)),
if member (false, map ('numberp,float (rest (pqxL)))) then
mydefint2_symb_lim (pqe,pqxL,pqyL)
else mydefint2_numer_lim (pqe,pqxL,pqyL))$
display2d:false$
ratprint:false$
---------------------------------------
Ted