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