Possible way to plot the solution density of diophantine equations

Question

Well, I'm trying to investigate the density of solutions to Diophantine equations. What is a general method to describe that density?

I've two functions:

$$\varphi\left(\text{a},\text{b},\text{c}\right)=\theta\left(\text{a},\text{b},\text{c}\right)\tag1$$

I will look for solutions in the following range: $\text{a}=\left\{\text{a}_0,\dots,\text{a}_\text{n}\right\}$ and $y=\left\{\text{b}_0,\dots,\text{b}_\text{m}\right\}$. So I will choose $\text{a}_0$ as starting value and let $\text{b}$ run from $\text{b}_0$ until $\text{b}_\text{m}$ and check if there follows a $\text{c}\in\mathbb{Z}$. Now, if I let $\text{b}_0$ and $\text{b}_\text{m}$ constant what should I plot on the x-axis if I want to plot the solution density on the y-axis?

---

Example and my work:

I've the Diophantine equation:

$$x^2+\left(\frac{y}{3}\right)^2=z\tag2$$

In order to solve it $x$, $y$ and $z$ has to be integers.

I will look for solutions in the following range: $x=\left\{5,\dots,7\right\}$ and $y=\left\{-10,\dots,10\right\}$. Now I found $\text{p}=21$ soltuions.

I used Mathematica to check them:

So, I should say that the density is given by:

$$\rho=\frac{\text{# solutions}}{\text{total possible solutions}}=\frac{21}{\left(1+7-5\right)\left(1+10-\left(-10\right)\right)}=\frac{21}{63}=\frac{1}{3}\tag3$$

Question: the general way of finding the solution density can be written as:

$$\rho=\frac{\text{# solutions}}{\left(1+x_\text{n}-x_0\right)\left(1+y_\text{m}-y_0\right)}\tag4$$

If I want to plot the function of $\rho$ on the y-axis what should be a smart choice to put on the x-axis? Assuming that I let $y_0$ and $y_\text{m}$ constant.

Answer 1

Clear["Global`*"]

ρ[x0_Integer, xn_Integer, y0_Integer, yn_Integer] :=
 Module[{x, y, z},
  Length[Solve[{x^2 + (y/3)^2 == z, x0 <= x <= xn, y0 <= y <= yn}, 
  {x, y, z}, Integers]]/((1 + xn - x0) (1 + yn - y0))]

Using symmetric bounds for x and y

DiscretePlot3D[ρ[-x, x, -y, y], {x, 1, 20}, {y, 1, 20}, 
 AxesLabel -> (Style[#, 14, Bold] & /@ {"x", "y", "\nρ "}),
 ColorFunction -> 
  Function[{x, y, ρ}, Piecewise[{{Red, y == 10}}, Blue]],
 ColorFunctionScaling -> False]

Fixing y to the interval {-10, 10} (colored Red above)

DiscretePlot[ρ[-x, x, -10, 10], {x, 1, 20},
 AxesLabel -> (Style[#, 14, Bold] & /@ {"x", "ρ"})]

Eliminating the negative value for x roughly cuts the number of solutions in half; however, the possible solution space is corresponding reduced by the same amount.

DiscretePlot[ρ[0, x, -10, 10], {x, 1, 20},
 AxesLabel -> (Style[#, 14, Bold] & /@ {"x", "ρ"})]