Four, (20 points) Given positive real numbers $a$ and $b$, it is known that real numbers $x$ and $y$ satisfy $a x^{2} - b x y + a y^{2} = 1$. Try to find the range of the bivariate function $f(x, y) = x^{2} + y^{2}$.

Translate the above text into English, please retain the line breaks and format of the source text, and output the translation result directly.
Let $x=\frac{m+n}{\sqrt{a}}, y=\frac{m-n}{\sqrt{a}}$. Then $a \cdot \frac{(m+n)^{2}}{a}-b \cdot \frac{m^{2}-n^{2}}{a}+a \cdot \frac{(m-n)^{2}}{a}=1$,

which means
$(2 a-b) m^{2}+(2 a+b) n^{2}=a$.
Thus, $f(x, y)=x^{2}+y^{2}$
$=\frac{(m+n)^{2}+(m-n)^{2}}{a}=\frac{2}{a}\left(m^{2}+n^{2}\right)$.
When $2 a-b>0$, i.e., $b<2 a$, from equation (1) we get
$(2 a+b)\left(m^{2}+n^{2}\right)=a+2 b m^{2} \geqslant a$,
$(2 a-b)\left(m^{2}+n^{2}\right)=a-2 b n^{2} \leqslant a$.
Therefore, $\frac{a}{2 a+b} \leqslant m^{2}+n^{2} \leqslant \frac{a}{2 a-b}$.
Thus, $\frac{2}{2 a+b} \leqslant \frac{2}{a}\left(m^{2}+n^{2}\right) \leqslant \frac{2}{2 a-b}$.
At this point, the range of $f(x, y)$ is
$\left[\frac{2}{2 a+b}, \frac{2}{2 a-b}\right]$.
When $2 a-b \leqslant 0$, i.e., $b \geqslant 2 a$, from equation (1) we get
$$
(2 a+b) n^{2}=a+(b-2 a) m^{2} \geqslant a \text {, }
$$

which means $n^{2} \geqslant \frac{a}{2 a+b}$.
Also, from $b-2 a \geqslant 0$ and equation (2), we know that $m^{2}$ can take any non-negative real number. Therefore,
$$
\begin{array}{l}
f(x, y)=\frac{2}{a}\left(m^{2}+n^{2}\right) \\
=\frac{2}{a}\left[m^{2}+\frac{a+(b-2 a) m^{2}}{2 a+b}\right] \\
=\frac{2}{2 a+b}+\frac{2}{a} \cdot \frac{2 b}{2 a+b} \cdot m^{2} \\
\geqslant \frac{2}{2 a+b} .
\end{array}
$$

At this point, the range of $f(x, y)$ is
$$
\left[\frac{2}{2 a+b},+\infty\right) \text {. }
$$

In summary, when $b<2 a$, the range of $f(x, y)=x^{2}+y^{2}$ is $\left[\frac{2}{2 a+b}, \frac{2}{2 a-b}\right]$; when $b \leqslant 2 a$, the range of $f(x, y)=x^{2}+y^{2}$ is $\left[\frac{2}{2 a+b},+\infty\right)$.