Functionality and Reliability II

*To*: mathgroup at smc.vnet.net*Subject*: [mg35495] Functionality and Reliability II*From*: Name <mee at home.com.redline.ru>*Date*: Tue, 16 Jul 2002 04:50:02 -0400 (EDT)*Sender*: owner-wri-mathgroup at wolfram.com

Actually there are cases where giving formal closed form answers for divergent sums/products leads to errors. As an illustration we can take two products, one diverging to infinity and another to zero: multiplying formal answers for Product[k(k+3)/(k+1)^2,{k,2,Infinity}] and Product[E^(-1/k),{k,2,Infinity}] gives an incorrect result for Product[k(k+3)/(k+1)^2 E^(-1/k),{k,2,Infinity}]. For sums/products that are done by matching against a pattern, it would be better to store the convergence conditions together with the pattern. Consider the example of taking integrals such as Integrate[1/(x^2-1),{x,-I,I}]. Problems of this type are doable with code such as given below: Clear[int] Options[int] = {PrincipalValue -> False, Assumptions -> {}}; int[ffz_, {z_Symbol, z1_, z2_}, opt___?OptionQ] := Module[ {init, main, jump, eq, ls, lse, sn, pvF, assum, ans}, SetAttributes[init, HoldAll]; init[fz_, pz_, qz_, Lroot_] := Module[ {Q}, fz = Cancel[Together[fz], Extension -> Automatic]; pz = Numerator[fz]; qz = Denominator[fz]; If[! PolynomialQ[pz, z] || ! PolynomialQ[qz, z], Return[$Failed] ]; fz = PolynomialQuotient[pz, qz, z]; pz = PolynomialRemainder[pz, qz, z]; Q = Evaluate[qz /. z -> #] &; Lroot = Array[Root[Q, #] &, Exponent[qz, z] ]; If[! MatchQ[MapThread[eq, {Lroot, RotateLeft[Lroot]} ], {(True | False) ...}], Return[$Failed] ]; Lroot = {#[[1]], Length[#]} & /@ Split[Sort[Lroot], eq]; ]; main[] := Module[ {fz = ffz, pz, qz, Lroot, Lroot2, ad, jmp, a0, z0, deg, tmp, A, i, j}, If[init[fz, pz, qz, Lroot] === $Failed, Return[$Failed] ]; ad = Integrate[fz, z]; jmp = 0; a0 = Coefficient[qz, z, Exponent[qz, z]]; For[i = 1, i <= Length[Lroot], i++, {z0, deg} = Lroot[[i]]; For[j = 1, j <= deg, j++, If[deg == 1, A = pz/D[qz, z] /. z -> z0, Lroot2 = Delete[Lroot, i]; tmp = Times @@ ((z - #[[1]])^#[[2]] & /@ Lroot2); A = D[pz/tmp, {z, deg - j}]/(a0(deg - j)!) /. z -> z0 ]; If[eq[A, 0], Continue[] ]; If[j == 1, ad += A Log[z - z0], ad += -A/((j - 1)(z - z0)^(j - 1))]; tmp = jump[z0, j]; If[tmp === $Failed, Return[$Failed] ]; jmp += A tmp ] ]; (ad /. z -> z2) - (ad /. z -> z1) - jmp ]; jump[z0_, deg_] := Module[ {x0, y0, x1, y1, x2, y2, x}, {x0, y0} = {Re[z0], Im[z0]}; {x1, y1} = {Re[z1], Im[z1]}; {x2, y2} = {Re[z2], Im[z2]}; If[! MatchQ[ {eq[x1, x0], eq[x2, x0], eq[y1, y0], eq[y2, y0]}, {(True | False) ..}], Return[$Failed] ]; If[eq[z1, z0] || eq[z2, z0], Return[$Failed] ]; If[deg != 1, If[eq[Abs[z0 - z1] + Abs[z2 - z0] - Abs[z2 - z1], 0], If[pvF && OddQ[deg], Return[0], Return[$Failed] ], Return[0], Return[$Failed] ] ]; If[eq[y1, y0] && eq[y2, y0], If[ls[x1, x0] && ls[x2, x0] || ls[x0, x1] && ls[x0, x2], Return[0] ]; If[pvF, Return[\[Pi] I sn[x1 - x2] ] ]; Return[$Failed] ]; If[ls[y1, y0] && ls[y2, y0] || lse[y0, y1] && lse[y0, y2], Return[0] ]; x = (y0 - y1)(x2 - x1)/(y2 - y1) + x1; If[ls[x0, x], Return[0] ]; If[ls[x, x0], Return[2\[Pi] I sn[y2 - y1] ] ]; If[eq[x, x0] && pvF, Return[\[Pi] I sn[y2 - y1] ] ]; $Failed ]; eq = Simplify[Chop[#1 - #2] == 0, assum] &; ls = Simplify[Chop[#1 - #2] < 0, assum] &; lse = Simplify[Chop[#1 - #2] <= 0, assum] &; sn = Simplify[Sign[Chop[#]], assum] &; {pvF, assum} = {PrincipalValue, Assumptions} /. Flatten[{opt, Options[int]}]; ans = main[]; ans /; ans =!= $Failed ] (*Examples: int[(x^5 - x^3 + 4*x)/(x^2 - 1)^3, {x, -1 - I, -1 + I}] int[(x^5 - x^3 + 4*x)/(x^2 - 1)^3, {x, -1 - I, -1 + I}, PrincipalValue -> True] // ComplexExpand int[(x^5 - x^3 + 4*x)/(x^2 - 1.`20)^3, {x, -1 - I, -1 + I}] int[(x^5 - x^3 + 4*x)/(x^2 - 1.`20)^3, {x, -1 - I, -1 + I}, PrincipalValue -> True] *) So far such or similar algorithm isn't implemented in Mathematica, and this isn't a problem per se. But the fact that Mathematica tries to handle such examples without having a robust algorithm is surprising. The situation with examples like Integrate[Min[Sin[x],2Sin[2x]],{x,0,Pi}] is similar: to deal with this problem we need to solve trigonometric equations with respect to the period, while Solve only finds some of the solutions. Of course the results of Integrate turn out to be incorrect. That doesn't mean having trigonometric solver is absolutely necessary. Things like a decent debugger or OpenGL support are more important. But here we are talking about fixing bugs, or rather about not creating bugs. In many other cases, as well as in the already discussed example of Limit[Sqrt[-1+a x],x->0] (if there are no means to give a conditional result or specify assumptions for a), wouldn't it be better just to leave the input unevaluated? As they say, 'don't pick at it, you'll only make it worse'. The key point is not that Mathematica can do better and cover wider range of problems, but rather that it should avoid giving results that only may happen to be correct. It seems that a large proportion of Mathematica's bugs should be attributed to that policy. Maxim Rytin m.r at prontomail.com