How Can I use Solve/Reduce Output

Question

Suppose I want x and y to be rationals

Solve[ x^2 + y^2 == 1, {x, y}, Rationals]

I am told:

Solve::svars: Equations may not give solutions for all "solve" variables. 

and given a "solution"

{{y -> ConditionalExpression[-Sqrt[ 1 - x^2], (x | y) ∈ Rationals && -1 <= x <= 1]},
 {y -> ConditionalExpression[ Sqrt[ 1 - x^2], (x | y) ∈ Rationals && -1 <= x <= 1]}}

Try Reduce

Reduce[ x^2 + y^2 == 1, {x, y}, Rationals]

Mathematica returns

(x | y) ∈ Rationals && -1 <= x <= 1 && (y == -Sqrt[1 - x^2] || y == Sqrt[1 - x^2])

Not much more help. Next try to see what Mathematica can do make a list of solutions:

FindInstance[ x^2 + y^2 == 1, {x, y}, Rationals, 10]

I am told that:

FindInstance::fwsol: Warning: FindInstance found only 2 instance(s), but it was not able
to prove 10 instances do not exist.

and given the two solutions found

{{x -> -1, y -> 0}, {x -> 1, y -> 0}}

Mathematica did not help too much here - but does show me that at least some solutions exist and that there may be more.

• Is this all I am going to get out of Mathematica in this case?
• Is this an answer that a mathematican wants to know and Mathematica has done its work?
• If I better understood Mathematica, could I use its other functions on the solution from Solve or Reduce to help me learn more about the solutions? (or do I just have to get my hands dirty and do the work myself - in this example it is pretty easy to find the general solution by hand)
• And maybe more generally, which Mathematica functions are able to do something more with this type of output from Solve or Reduce? (e.g. can I plot it?)

Answer 1

Direct solutions using Solve

Since there are infinitely many solutions to the problem at hand we prefer to choose a bound on the solution space, i.e. we solve (x/y)^2 + (z/w)^2 == 1 in Integers and select only solutions satisfying x/y <= z/w and 0 < x <= y <= 50 and z <= w <= 50. It takes a bit but we know we have found all solutions satisying given criteria.

Union[{#1/#2, #3/#4} & @@@ Solve[(x/y)^2 + (z/w)^2 == 1 && 0 < x <= y <= 50 &&
                                  x/y <= z/w && z <= w <= 50, {x, y, z, w},
                                                                  Integers][[All, All, 2]]]

{{9/41, 40/41}, {7/25, 24/25}, {12/37, 35/37}, {5/13, 12/13}, {8/17, 15/17},
 {3/5, 4/5}, {20/29, 21/29}}

Now we can check if these are solutions :

And @@ (#1^2 + #2^2 == 1 & @@@ %)

 True

Well, indeed.

Edit 1

Solutions generated by Table of Pythagorean triples

Some more efficient way of generating rational solutions of the equation x^2 + y^2 == 1 is based on Pythagorean triples. We can write a table of numbers of the form e.g. $(\frac{x^2-y^2}{x^2+y^2},\frac{2 x y}{x^2+y^2})$ and select only distinct ones :

pts = Union @ Flatten[Table[{(x^2 - y^2)/(x^2 + y^2), (2 x y)/(x^2 + y^2)},
                            {y, 20}, {x, 20}], 1];
Length @ pts

255    

ParametricPlot[ {Cos[u], Sin[u]}, {u, 0, Pi/2}, PlotStyle -> Thick, 
                                   Epilog -> { Red, PointSize[0.006], Point @ pts}]

Solutions generated by Fibbonacci numbers

We could find only special solutions making use of observations by Dujella that Fibbonacci numbers $F_{n}$ generate distinct Pythagorean triples $(F_{n} F_{n+3},2F_{n+1}F_{n+2},F_{n+1}^2+F_{n+2}^2)$ although not exhaustively. Thus proceeding this way we can find arbitrary many distinct solutions quite effectively. We find the first 100 such solutions :

ptsF = 
  Table[{( Fibonacci[n] Fibonacci[n + 3])/( Fibonacci[n + 1]^2 + Fibonacci[n + 2]^2),
         ( 2Fibonacci[n + 1] Fibonacci[n + 2])/( Fibonacci[n + 1]^2 + Fibonacci[n + 2]^2)},
         {n, 100}];

And @@ (#1^2 + #2^2 == 1 & @@@ ptsF)

True

Let's write some of them :

ptsF[[;; 20]]

ParametricPlot[{Cos[u], Sin[u]}, {u, 0, Pi/2}, Epilog -> {Red, Point @ ptsF}]

Edit 2

Since there could be some kind of misunderstanding I add another way to visualize the solution set for various numbers of solutions constructed with tables of Pythagorean triples.

Clear[pt]
pt[k_Integer /; OddQ[k] && k > 0] := 
  Union @ Flatten[ 
            Table[{(x^2 - y^2)/(x^2 + y^2), (2 x y)/(x^2 + y^2)}, 
                                       {y, -k, k, 2}, {x, -k, k} ], 1]

Now we check what is distribution of solution points in various cases, e.g.

GraphicsRow[ Table[ 
                    ParametricPlot[ {Cos[u], Sin[u]}, {u, 0, 2 Pi},
                                                      PlotStyle -> Thin, 
                                                      Epilog -> { Red, Point @ p} ],
                   {p, {pt[9], pt[17], pt[55]}}]]