Menu
▾
▴
[Maxima-commits] CVS: maxima/share/contrib/integration abs_integrate.html, 1.7, 1.8 abs_integrate.mac, 1.29, 1.30 abs_integrate.texi, 1.8, 1.9 hyperint.mac, 1.4, 1.5 rtest_abs_integrate.mac, 1.29, 1.30
|
From: Barton W. <wil...@us...> - 2011-01-13 17:46:59
|
Update of /cvsroot/maxima/maxima/share/contrib/integration In directory sfp-cvsdas-4.v30.ch3.sourceforge.com:/tmp/cvs-serv27208 Modified Files: abs_integrate.html abs_integrate.mac abs_integrate.texi hyperint.mac rtest_abs_integrate.mac Log Message: o Fix for SF bug abs_integrate of unit step - ID: 3153169 o Fix for SF bug abs_integrate --> errors from apply - ID: 3154351 o Fix for abs_integrate + composisiton of unit step - ID: 3153170 o Append regression tests for bugs 3154351 & 3153170 o Update user docummentation (removed make_dummy from user documentation) o Change make_dummy to new_variable (changes for hyperint.mac) o Update regression tests (signum is now fully multiplicative). Index: abs_integrate.html =================================================================== RCS file: /cvsroot/maxima/maxima/share/contrib/integration/abs_integrate.html,v retrieving revision 1.7 retrieving revision 1.8 diff -u -d -r1.7 -r1.8 --- abs_integrate.html 11 Aug 2010 19:06:30 -0000 1.7 +++ abs_integrate.html 13 Jan 2011 17:46:49 -0000 1.8 @@ -1,6 +1,6 @@ <HTML> <HEAD> -<!-- Created by texi2html 1.56k from abs_integrate.texi on 11 August 2010 --> +<!-- Created by texi2html 1.56k from abs_integrate.texi on 8 January 2011 --> <TITLE>Untitled Document</TITLE> </HEAD> @@ -460,32 +460,8 @@ </DL> -<P> -<DL> -<DT><U>Function:</U> <B>make_dummy(e,x)</B> -<DD><A NAME="IDX9"></A> - - -<P> -Append "%" to the symbol <EM>x</EM> until <EM>x</EM> is <I>not</I> a variable in the expression <EM>e</EM>. - - - -<PRE> - (%i1) make_dummy((x-y)*x, x); - (%o1) x% - - (%i2) make_dummy((x-y)*x%, x); - (%o2) x%% -</PRE> - -<P> -<B>To use</B> <TT>`load(abs_integrate)'</TT> - - -</DL> <P><HR><P> -This document was generated on 11 August 2010 using +This document was generated on 8 January 2011 using <A HREF="http://wwwinfo.cern.ch/dis/texi2html/">texi2html</A> 1.56k. </BODY> </HTML> Index: abs_integrate.mac =================================================================== RCS file: /cvsroot/maxima/maxima/share/contrib/integration/abs_integrate.mac,v retrieving revision 1.29 retrieving revision 1.30 diff -u -d -r1.29 -r1.30 --- abs_integrate.mac 11 Aug 2010 19:06:30 -0000 1.29 +++ abs_integrate.mac 13 Jan 2011 17:46:49 -0000 1.30 @@ -1,5 +1,5 @@ /* - Copyright 2008, 2009, 2010 by Barton Willis + Copyright 2008, 2009, 2010, 2011 by Barton Willis Maxima code for integration of signum, abs, max, min, and unit_step. This is free software; you can redistribute it and/or @@ -34,9 +34,9 @@ while l # [ ] do ( i : first(l)(q,x), l : rest(l), - if (i # false) then ( + if (i # false and freeof('integrate, nounify('integrate), i)) then ( ii : i, /* accept i as a simplification of ii. */ - if freeof('integrate, i) then l : [])), + l : [])), ii); extra_definite_integration_methods : ['abs_defint]; @@ -63,15 +63,21 @@ /* Convert abs, max, min, and unit_step to signum expressions */ convert_to_signum(e) := block([opsubst : true], e : unit_step_mult_simp(e), - e : subst('unit_step = lambda([s], (signum(s) + 1) / 2), e), + e : subst('unit_step = lambda([s], (signum(s)*(signum(s)+1))/2), e), e : subst('max = lambda([[x]], xreduce(lambda([a,b], (a + b +abs(a-b))/2),x)),e), e : subst('min = lambda([[x]], xreduce(lambda([a,b], (a + b -abs(a-b))/2),x)),e), subst('abs = lambda([s], s * signum(s)), e)); /* Do signum(a) * signum(b) --> signum(a * b). Maybe it should, but when n is a positive integer, - apply_signum_mult doesn't do signum(a)^n --> signum(a^n).*/ - + apply_signum_mult doesn't do signum(a)^n --> signum(a^n). + Starting with Maxima 5.23, the general simplifier does signum(a * b) --> signum(a) * signum(b). + Often (but not always!) the general simplifier reverts what this function does. Why not + always? We have apply_signum_mult_id(signum(x) * signum(1/x)) --> 1, but the general simplifier + doesn't change signum(x) * signum(1/x). + + This function isn't used--maybe it will eventually be removed. */ + apply_signum_mult_id (e) := subst(["*" = lambda([[x]], block([p : [], q : []], for xk in x do ( @@ -92,21 +98,13 @@ if freeof('max, 'min, 'inf, p) then q * unit_step(p) else e) else map('unit_step_mult_simp, e)); -/* Do signum(a) * signum(b) --> signum(a * b) and e * signum(e) --> abs(e) */ +/* Do e * signum(e) --> abs(e) */ signum_to_abs(e) := block([l], - e : apply_signum_mult_id(e), l : flatten(gatherargs(e,'signum)), for lk in l do ( e : ratsubst(abs(lk), lk * signum(lk),e)), e); -/* Do signum(a*b) --> signum(a) * signum(b) */ - -signum_expand(e) := block([sublis_apply_lambda : true], - sublis(['signum = lambda([s], - (s : factor(s), - if safe_op(s) = "*" then xreduce("*", map('signum, args(s))) else signum(s)))], e)); - /* If e has the form w * signum(p1(x)) * signum(p2(x)) * ... * signum(pn(x)), where w is either freeof signum terms or free of x, return [w, append(q1,q2,...,qn)], where qk are the factors of pk; otherwise, return false. */ @@ -129,6 +127,7 @@ signum_int(q,x) := block([w : 1, acc : [], sgn, v, f, listconstvars : true], q : convert_to_signum(q), + q : almost_everywhere_simp(q), v : listofvars(q), q : block([expandwrt_denom : true], expandwrt(q,'signum)), if (f : abs_int_extra(q,x)) # false then f @@ -141,10 +140,13 @@ signum_int_helper(signum(sgn) * w, sort(acc), x))) else if (safe_op(q) = "+") then ( f : map(lambda([s], signum_int(s,x)), args(q)), - if every(lambda([s], s # false), f) then xreduce("+", f) else false) + if member(false, f) then false else xreduce("+", f)) else block([extra_integration_methods : [], extra_definite_integration_methods : []], integrate(q,x))); +/* This function is deprecated. Either use new_variable (from to_poly_solve_extra) or +gensym. */ + make_dummy(e,x) := block([listconstvars : true, v : listofvars(e)], while member(x, v) do x : concat(x,"%"), x); @@ -155,7 +157,7 @@ if emptyp(rest(l)) then ( f : integrate(q,x), if not(freeof('integrate, f)) then ( - v : make_dummy(q,x%), + v : new_variable('general), f : funmake(nounify('integrate), [subst(x = v, q), v, first(l), x])), signum(x - first(l)) * (f - subst(x = first(l), f))) else ( @@ -174,20 +176,27 @@ xo : xi), acc); -partitions_interval_p(l,[endpts]) := ( - if listp(l) then ( +partitions_interval_p(l,[endpts]) := block([ok,lo,hi], + ok : if length(endpts) = 2 then ( + [lo,hi] : endpts, + every(lambda([s], (is(lo <= s) = true) and (is(s <= hi) = true)), l)) + else endpts = [], + if listp(l) and ok then ( l : listify(setify(append(endpts, l))), - l : sort(l, lambda([a,b], is(a < b) = true)), + l : sort(l, lambda([a,b], csign(b-a) = 'pos)), if every(lambda([a,b], is(a <= b)=true), l, endcons('inf, rest(l))) then l else false) else false); -abs_defint(e,x,lo,hi) := block([f,l], +abs_defint(e,x,lo,hi) := block([f,l, prederror : false], if is(hi < lo)=true then -abs_defint(e,x,hi,lo) else if is(lo = hi) then 0 else ( + e : simp_assuming(e,lo <= x, x <= hi), e : convert_to_signum(e), + e : almost_everywhere_simp(e), l : xreduce('append, (gatherargs(e, 'signum))), l : flatten(map(lambda([s], real_linear_factors(s,x)), l)), /* flatten([[1], false]) --> [1,false] */ + l : sublist(l, lambda([s], lo <= s and s <= hi)), l : delete(false,l), l : partitions_interval_p(l,lo,hi), if l # false then ( @@ -200,6 +209,20 @@ simp_assuming(e, [fcts]) ::= buildq([e,fcts], unwind_protect((apply(assume, fcts), expand(e,0,0)), apply(forget,fcts))); +/* Replace signum(z)^2 --> 1, where z is nonzero and a polynomial. For a nonzero polynomial z, +signum(z)^2 -1 is zero everywhere except at each zero of z. The almost_everywhere_simp simplification +is OK for an integrand--thus integrate(e,x) = integrate(almost_everywhere_simp(e),x). + +This function could be extended to do signum(z) --> 1, when z is a nonzero polynomial at csign(z) = 'pn, +for example. Another extension is %if(x=0,xxx,zzz) --> zzz. I'll save these for another day.*/ + +almost_everywhere_simp(e) := block([l], + l : sublist(flatten(gatherargs(e,'signum)), lambda([s], s # 0 and polynomialp(s,listofvars(s)))), + e : ratexpand(e), + for lk in l do ( + e : ratsubst(1, signum(lk)^2,e)), + e); + /* For integrands of the form F(x, |x-c|), integrate F(x,-x+c) and F(x,x-c). Return a signum expression that is continuous at c; when the integrand doesn't have this form, return false. The error catch on integrate is needed: try integrate(1/(x + abs(x)),x), for example. @@ -239,7 +262,7 @@ k : flatten(gatherargs(q,'signum)), if k # [] and every(lambda([s], polynomialp(s,[x], 'numberp, lambda([k],is(k=0 or k=1)))), k) then ( k : map(lambda([s], rhs(first(linsolve(s,x)))),k), - %x : make_dummy(q,x), + %x : new_variable('general), q : sublis([x = %x],q), k : listify(setify(k)), k : partitions_interval_p(k), @@ -277,8 +300,8 @@ l : union(l, setify(xreduce('append, gatherargs(q,'ceiling)))), l : subset(l, lambda([s], not(freeof(x, s)))), if l = set(x) then ( - fl : make_dummy([q,x], 'fl), - cl : make_dummy([q, x, fl], 'cl), + fl : new_variable('general), + cl : new_variable('general), q : subst([floor(x) = fl, ceiling(x) = cl], q), q : integrate(q, x), if freeof('integrate, q) then ( Index: abs_integrate.texi =================================================================== RCS file: /cvsroot/maxima/maxima/share/contrib/integration/abs_integrate.texi,v retrieving revision 1.8 retrieving revision 1.9 diff -u -d -r1.8 -r1.9 --- abs_integrate.texi 11 Aug 2010 19:06:30 -0000 1.8 +++ abs_integrate.texi 13 Jan 2011 17:46:49 -0000 1.9 @@ -322,20 +322,3 @@ @end deffn -@deffn {Function} make_dummy(e,x) - -Append "%" to the symbol @math{x} until @math{x} is @i{not} a variable in the expression @math{e}. - -@example - (%i1) make_dummy((x-y)*x, x); - (%o1) x% - - (%i2) make_dummy((x-y)*x%, x); - (%o2) x%% -@end example - -@b{To use} @file{load(abs_integrate)} - - - -@end deffn \ No newline at end of file Index: hyperint.mac =================================================================== RCS file: /cvsroot/maxima/maxima/share/contrib/integration/hyperint.mac,v retrieving revision 1.4 retrieving revision 1.5 diff -u -d -r1.4 -r1.5 --- hyperint.mac 21 Jan 2010 11:36:15 -0000 1.4 +++ hyperint.mac 13 Jan 2011 17:46:49 -0000 1.5 @@ -67,7 +67,7 @@ /* Express sigma as a linear combination of the members of kern. */ xresidue(sigma, kern, x) := block([zip, %g, i, n, k], - %g : make_dummy([sigma, kern], '%g), + %g : new_variable('general), n : length(kern), zip : partfrac(sigma - kern . makelist(concat(%g,i),i,1,n),x), zip : if safe_op(zip) = "+" then args(zip) else [zip], @@ -92,8 +92,8 @@ if kern # false then ( kern : listify(kern), sigma : partfrac(ratsimp(diff(e,x) / e),x), - a : make_dummy([e,x],'a), - b : make_dummy([e,x],'b), + a : new_variable('general), + b : new_variable('general), dkern : map(lambda([s], diff(s,x) / s), kern), r : xresidue(sigma, dkern, x), @@ -114,7 +114,7 @@ f : block([?errorsw : true], errcatch(subst(sk, f))), if f # [] then ( f : first(f), - %fo : make_dummy([f,e,x], '%fo), + %fo : new_variable('general), f : %fo * f, cnst : subst(sk, diff(f,x) * f^a * (1 - f)^b / e), cnst : infapply(cnst,lambda([s], radcan(rootscontract(s)))), @@ -164,7 +164,7 @@ if kern # false then ( kern : listify(kern), sigma : partfrac(ratsimp(diff(e,x) / e),x), - %k : make_dummy([e,x],'%k), + %k : new_variable('general), dkern : map(lambda([s], diff(s,x) / s), kern), r : xresidue(sigma, dkern, x), l : if r = [ ] then [] else generate_elliptic_candidates(dkern, r, %k), @@ -193,7 +193,7 @@ f : block([?errorsw : true], errcatch(subst(sk, f))), if f # [] then ( f : first(f), - %fo : make_dummy([f,e,x], '%fo), + %fo : new_variable('general), f : %fo * f, cnst : subst(sk, (diff(f,x) / (sqrt(1-f^2) * (sqrt(1-%k^2 * f^2)))) / e), cnst : infapply(cnst,lambda([s], radcan(rootscontract(s)))), Index: rtest_abs_integrate.mac =================================================================== RCS file: /cvsroot/maxima/maxima/share/contrib/integration/rtest_abs_integrate.mac,v retrieving revision 1.29 retrieving revision 1.30 diff -u -d -r1.29 -r1.30 --- rtest_abs_integrate.mac 6 Oct 2010 19:21:43 -0000 1.29 +++ rtest_abs_integrate.mac 13 Jan 2011 17:46:49 -0000 1.30 @@ -42,8 +42,8 @@ i : integrate(abs(x) + abs(1-x),x); x*abs(x)/2+(x-1)*abs(x-1)/2$ -factor(signum_to_abs(diff(%,x))); -abs(x)+abs(x-1)$ +factor(diff(i,x) - (abs(x) + abs(1-x))); +0$ (remvalue(i),0); 0$ @@ -70,13 +70,10 @@ (x^2*signum(x))/4+(3*x^2)/4$ i : integrate(max(x,x^2),x); -(((2*x^3-3*x^2)/12+1/12)*signum(x-1)+1/12)*signum(x)+(x^3/3+x^2/2)/2$ - -di : simp_assuming(diff(i,x), notequal(x,0),notequal(x,1)); -(6*x^2-6*x)*signum(x-1)*signum(x)/12+(x^2+x)/2$ +(((2*x^3-3*x^2)/12+1/12)*signum(x-1)+1/12)*signum(x)+x^3/6+x^2/4$ -makelist(ratsimp(subst(x = k, max(x,x^2) - di)),k, [-1066, -1, 0, 1, %pi, sqrt(2), 3/2, 42, 1932]); -[0,0,0,0,0,0,0,0,0]$ +factor(convert_to_signum(diff(i,x) - max(x,x^2))); +0$ limit(i,x,1, 'minus) - limit(i,x,1, 'plus); 0$ @@ -108,8 +105,8 @@ integrate(unit_step(x^2-1),x); ((x/2-1/2)*signum(x-1)-1)*signum(x+1)+x/2$ -expand(signum_expand(convert_to_signum(diff(%,x) - unit_step(x^2-1)))); -0$ +almost_everywhere_simp(scanmap('factor, convert_to_signum(diff(%,x) - unit_step(x^2-1)))); +0$ integrate(abs(x),x,1,2); 3/2$ @@ -135,10 +132,11 @@ signum_to_abs(integrate(unit_step(x),x,a,b)); (abs(b)+b)/2-(abs(a)+a)/2$ -(i : integrate(f(x) * abs(x-a),x), ev(i, f(x) := x, 'integrate)); -((2*x^3-3*a*x^2)/6+a^3/6)*signum(x-a)$ +/* With hyper_int, we get an interesting result... */ +block([extra_integration_methods : ['signum_int]], i : integrate(f(x) * abs(x-a),x), ev(i, f(x) := x, 'integrate)); +x^2*abs(x-a)/2-(2*x^3/3-2*a^3/3)*signum(x-a)/4$ -signum_to_abs(simp_assuming(diff(%,x),notequal(x,a))); +factor(simp_assuming(diff(signum_to_abs(%),x),notequal(x,a))); x*abs(x-a)$ (i : integrate(f(x) * max(x,-x),x,0,1), ev(i, f(x) := x,integrate)); @@ -147,11 +145,11 @@ subst([a=-1,b=2/3, c=1], integrate(abs((x-a)*(x-b)*(x-c)),x,b,c)); 11/972$ -(i : integrate(abs(x^2-1)*x,x),0); -0$ +i : integrate(abs(x^2-1)*x,x); +(x^2-1)*abs(x^2-1)/4$ -signum_to_abs(diff(i,x)); -x*abs(x^2-1)$ +factor(diff(i,x) - abs(x^2-1)*x); +0$ is(equal(limit(i,x,-1,'minus), limit(i,x,-1,'plus))); true$ @@ -186,11 +184,11 @@ integrate(abs(x-1) / (x + abs(x-5)),x,9,-7) + integrate(abs(x-1) / (x + abs(x-5)),x,-7,9); 0$ -integrate(signum(x) * f(x),x); -signum(x)*integrate(f(x%),x%,0,x)$ +nicedummies(integrate(signum(x) * f(x),x)); +'integrate(signum(x) * f(x),x)$ -integrate(signum(x-a) * f(x),x); -signum(x-a)*integrate(f(x%),x%,a,x)$ +nicedummies(integrate(signum(x-a) * f(x),x)); +'integrate(signum(x-a) * f(x),x)$ is(equal(op(integrate(signum(x + sqrt(1-x)),x)), nounify('integrate))); true$ @@ -406,13 +404,13 @@ 0$ integrate(1/ (1 + max(x^2,x)),x,-1,1); -log(2)+%pi/4$ +2*(log(4)/2-log(2)/2)+%pi/4$ integrate(x/ (1 + max(x^2,1)^2),x,minf,inf); 0$ integrate(x/ (1 + max(x^2,1)^2),x,0,inf); -%pi/8$ +%pi/8 + 1/4$ /* the call to expand is a workaround for an integrate bug */ expand(integrate(mod(x,1),x),0,0); @@ -432,7 +430,7 @@ /* SF bug ID: 3034415; the antidervative is non-ideal, but at least it nolonger signals an error. */ integrate(1/abs(x),x); -'integrate(1/(x*signum(x)),x)$ +'integrate(1/abs(x),x)$ (extra_integration_methods : cons('conditional_integrate,extra_integration_methods),0); 0$ @@ -483,7 +481,22 @@ 'integrate(abs(1 - a*z),z); %if(a#0,((a*z-1)*abs(a*z-1))/(2*a),z)$ -(print("time = ", elapsed_real_time () - start), is(elapsed_real_time () - start < 45)); +integrate(unit_step(x),x,99,100); +1$ + +abs_defint(unit_step(x),x,99,100); +1$ + +abs_defint(unit_step(unit_step(x)),x,-1,0); +0$ + +integrate(unit_step(unit_step(x)),x,-1,0); +0$ + +integrate(signum(abs(x)),x,-2,2), prederror : false; +4$ + +(print("time = ", elapsed_real_time () - start), is(elapsed_real_time () - start < 70)); true$ (remvalue(e,i,di,f1,start), remove(x,noninteger), remfunction(convolution, unit_box),0); |