Message Boards Message Boards

Why does FullSimplify not recognize this function evaluates to zero?

Hi all, In trying to compare a student solution to the "book" solution for a function, I want to use FullSimplify[] to form a difference function to see whether the two forms are identical. But for these two functions of x>0 and y > 0,

fbook[x_, y_] := Sqrt[((1 + x) (x + y))/y]
fstudent[x_, y_] := Sqrt[((1 + x) y)/(x + y)] + x Sqrt[(1 + x)/(x y + y^2)]
diff[x_, y_] := fstudent[x, y] - fbook[x, y]

I find that FullSimplify[diff[x,y], {x>0, y>0}] does not return zero... unless I multiply diff[x,y] either by Sqrt[y] or Sqrt[x+y].

I did read the documentation on Simplify and FullSimplify, and understand that both of these try to find a simpler form of their argument, with no guarantees. And I learned that sometimes even the order of variables affects the results. I am just curious if anyone can provide some insight into this case.

Also, how could/should one quickly verify that two functions are in fact one and the same (in terms of returning the same values for all valid inputs), given that FullSimplify seems to have trouble with this fairly straightforward example? Notebook attached.

Attachments:
POSTED BY: Philip Blanco
7 Replies

You're welcome. Let me know if you have any questions. To inspect the code, consider changing the Module[{...}, line to

Module[{},

with no localized variables. That way, after a test run you can look at the results in each variable. Also change npts = 10 or some other number, so it's comfortable to check testpts.

You might replace

testpts = Transpose[
     RandomReal[#, npts, WorkingPrecision -> precision] & /@ domains
     ];
   AllTrue[testpts, e1 == e2 /. Thread[vars -> #] &]

with

testpts = Transpose[
     Rationalize[RandomReal[#, npts], Max[Abs@#] * 10^-5] & /@ domains
     ];
   AllTrue[testpts, PossibleZeroQ[e1 - e2 /. Thread[vars -> #] &]

Rationalize is used because you need to pass exact values to e1 and e2 for PossibleZeroQ to work. The factor 10^-5 in the second argument (indirectly) limits the size of the denominator (and therefore the numerator, too) in the fractions produced. Big numbers in the fraction can lead to very long computations sometimes. I don't know if 10^-5 is the best choice, but it doesn't seem like it's probably a bad choice.

I suggest it because I think PossibleZeroQ is a little more rigorous than the simple equality check in my first version.

POSTED BY: Michael Rogers

Wow thank you MIchael for the detailed solution! I shall adapt and make use of that code (it will take me a little while to understand all the bells and whistles), since it looks like you have produced a very useful tool for teachers, and for other professionals who want to check different forms of functions that should produce the same answer over a domain. Much more reliable - and quantifiable - than plotting and eyeballing.

Oh and luckily for my field in education, "ill-conditioned problem" = "learning opportunity".

POSTED BY: Philip Blanco

Numerically checking equality is a practical way to be confident (to a degree) of a right/wrong answer. It is usually faster than symbolic analysis. There is a small probability that the check will fail. And numerical analysis comes into play. There are some ill-conditioned problems for which there is a high probability that a correct answer would be rejected. Mathematica has some sophisticated tools you can use if you wish the check to be virtually always correct. However, ill-conditioned problems might not come up very often in homework, and you might be willing to deal with computer mistakes on a case-by-case basis. This sort of checking needs to be viewed as a practical matter, and you need to decide what is worth worrying about and what isn't.

Simple numerical check

Here is a simple numerical check to test if two expressions are equal. Usage: simpleNEqual[e1, e2, {x, a, b}] or simpleNEqual[e1, e2, {x, a, b}, {y, c, d}] etc., with however many standard iterators {x, a, b} like in Integrate that you need, one for each variable. (This could be extended to regions, using RandomPoint instead of RandomReal, if rectangular domains won't work.)

Weaknesses:

  • It checks numerical equality at 100 points (npts = 100). It is possible that two unequal expressions agree on 100 random points. Possible but unlikely.
  • It uses arbitrary-precision numbers of 32 digits precision (precision = 32). It is possible that with round-off error, an expression might lose all digits of precision and the calculated value become unreliable. This happens with single-digit floating point number (~8 digits of precision) with enough frequency that when double-precision (~16 digits) became common in CPUs, the people rejoiced. However, everyone who does numerical computation finds occasionally that 16 digits is not enough. With 32 digits, it should be rare that precision loss leads to an unreliable result, but it can happen.

Code:

ClearAll[simpleNEqual];
simpleNEqual[e1_, e2_, dom : {_, _, _} ..] :=
  Module[{vars, domains, testpts, npts, precision},
   (* set up *)
   vars = {dom}[[All, 1]];
   domains = {dom}[[All, 2 ;; 3]] /. {-Infinity -> -1000, Infinity -> 1000}; (* impose finite limits on Infinity *)
   (* internal parameters - could be arguments/options *)
   npts = 100;     (* enough points? *)
   precision = 32; (* too low, might fail; too high, waste time *)
   (* actual checks *)
   testpts = Transpose[
     RandomReal[#, npts, WorkingPrecision -> precision] & /@ domains
     ];
   AllTrue[testpts, e1 == e2 /. Thread[vars -> #] &]
   ];

Examples:

simpleNEqual[fstudent[x, y], fbook[x, y], {x, 0, Infinity}, {y, 0, Infinity}] // AbsoluteTiming
(*  {0.0036, True}  *)

simpleNEqual[fstudent[x, y], (1 + 10^-29) fbook[x, y],   (* the small difference is detected *)
 {x, 0, Infinity}, {y, 0, Infinity}]
simpleNEqual[fstudent[x, y], (1 + 10^-30) fbook[x, y],   (* the difference is too small and is undetected *)
 {x, 0, Infinity}, {y, 0, Infinity}]
(*
  False
  True
*)

Sprucing the code up a bit. You could add options to control the value of npts and precision. You could also put in checks to see if there were problems and report via Sow the problem and the input that caused the problem; use Reap to get the information. For instance, define an error message and change the AllTrue[..] code as follows:

simpleNEqual::prec = "Less than two-digit precision at ``.";

AllTrue[testpts,
 With[{msg = $MessageList},
   Check[res = {e1, e2} /. Thread[vars -> #], (* check if evaluation of the expressions gives errors/warnings *)
    Sow[# -> Complement[$MessageList, msg], "error"]];
   If[Precision[res] < 2 && res[[1]] != 0 && res[[2]] != 0,  (* check if too much precision loss *)
    Message[simpleNEqual::prec, Sow[Thread[vars -> #], "precision"]]];
   res = Equal @@ res;
   If[! res, Sow[Thread[vars -> #], "unequal"]];             (* Sow the input that proves e1 != e2 *)
   res] &]
POSTED BY: Michael Rogers

Thank you for your detailed responses - I see that I have much to learn! I look forward to additional replies and suggestions. I am curious how teachers could check answers given in different forms in the general case.

(I think I heard that some online homework systems just evaluate the difference function for a range of possible valid arguments to check that it stays below some tolerance.)

POSTED BY: Philip Blanco

Simplify does not handle this example very well. However, you can try Reduce:

Reduce[diff[x, y] == 0, Reals]

and plot the region of validity:

RegionPlot[
 Evaluate@Reduce[diff[x, y] == 0, Reals], {x, -5, 5}, {y, -5, 5}]

(it does not show the 1-dimensional pieces). The diff[x,y] is zero only on a subset of its domain:

RegionPlot[
 Evaluate@Reduce@FunctionDomain[diff[x, y], {x, y}, Reals], {x, -5, 
  5}, {y, -5, 5}]

For example,

In[11]:= diff[-4, 2]

Out[11]= -2 Sqrt[3]
POSTED BY: Gianluca Gorni

Philip Blanco specified the domain "x>0 and y > 0". Perhaps you overlooked that, or discounted it for some reason?

Reduce handles some forms of the problem:

Reduce[x > 0 && y > 0 && diff[x, y] == 0]
(*  y > 0 && x > 0  *)

Reduce[Implies[x > 0 && y > 0, diff[x, y] == 0], Reals]
(*  True  *)

Reduce[ForAll[{x, y}, x > 0 && y > 0, diff[x, y] == 0]]
(*  True  *)

Reduce[ForAll[{x, y}, x > 0 && y > 0, z == diff[x, y]], Reals]
(*  z == 0  *)

Simplify@Reduce[Implies[x > 0 && y > 0, z == diff[x, y]], Reals]
(*  z == 0 || x <= 0 || y <= 0  *)

The alternatives x <= 0 and y <= 0 in the last form can be eliminated with Simplify[%, x > 0 && y > 0]; the same can be done with the first form to reduce it to True.

The forms above where Reals are specified seem to need the domain restriction. It appears that in those forms, Reduce cannot rule out, or otherwise deal with, cases such as x + y == 0 or x + y < 0, etc.

POSTED BY: Michael Rogers

1) The expanded radical by itself is not simpler:

    Simplify`SimplifyCount[Sqrt[((1 + x) (x + y))/y]]
    Simplify`SimplifyCount[Sqrt[1 + x] Sqrt[1/y] Sqrt[x + y]]
    (*
      15
      23
    *)

So if such a transformation is tried, it might be rejected out of hand or before it led to a simpler result. And it's possible that such a transformation won't be tried, because it is known that it doesn't simplify the term, even if further down the line it would cancel with something else. The side-alleys that can be pursued in simplification may be infinite and need to be limited in some way.

2) The transformation you seek is PowerExpand. However, because of the issue in 1) above, adding it does not solve your problem:

Simplify[diff[x, y], x > 0 && y > 0, TransformationFunctions -> {Automatic, PowerExpand}]
(*  Sqrt[((1 + x) y)/(x + y)] - Sqrt[((1 + x) (x + y))/y] + x Sqrt[(1 + x)/(x y + y^2)]  *)

3) The trick is to simplify after the expansion, before the PowerExpand transformation is rejected. Both of the following work:

Simplify[PowerExpand@diff[x, y], x > 0 && y > 0]
(*  0  *)

Simplify[diff[x, y], x > 0 && y > 0, TransformationFunctions -> {Automatic, Simplify@*PowerExpand}]
(*  0  *)

4) Here is a custom transformation for distributing the square roots, assuming they are not in the denominator. Note that Sqrt[u] expand to Power[u, 1/2] and 1/Sqrt[u] expands to Power[u, -1/2]. This is important when using pattern-matching, since the second will not match Sqrt[_].

sqrtExpand[e_] := e /. rad : Sqrt[Times[factors__]] :> 
    Thread[rad, Times] /;                             (* Thread[] distributes the Power[..,1/2] over Times *)
      Simplify[And @@ Thread[List @@ factors >= 0]];  (* The condition is that all factors are nonnegative *)
Assuming[x > 0 && y > 0,
 Simplify[diff[x, y], TransformationFunctions -> {Automatic, Simplify@*sqrtExpand}]
 ]
(*  0  *)

5) Note the using of Assuming[x > 0 && y > 0,..] instead of Simplify[.., x > 0 && y > 0]. Simplify does not add its assumptions to $Assumptions; instead it passes them to certain transformations it uses (apparently). However, it does use any assumptions in $Assumptions. Consequently, the simplification in Simplify@*sqrtExpand would not be done with the assumption x > 0 && y > 0 if I used Simplify[.., x > 0 && y > 0]; however, with Assuming, the assumptions will be used in all functions called during the computation that use $Assumptions.


Well, that is about all I can add. I realize it doesn't completely answer your question, and some of the remarks go beyond what you asked; but I hope it gives you some insight that you seek.

POSTED BY: Michael Rogers
Reply to this discussion
Community posts can be styled and formatted using the Markdown syntax.
Reply Preview
Attachments
Remove
or Discard

Group Abstract Group Abstract