The answer is $f(x) = C_1x + \frac{C_2}{x}$ with arbitrary constants $C_1$ and $C_2$.

Solution 1. Fix a real number $a > 1$, and take a new variable $t$. For the values $f(t), f(t^2), f(at)$ and $f(a^2t^2)$, the relation of the assumption provides a system of linear equations:
$$x = y = t : \left(t + \frac{1}{t}\right) f(t) = f(t^2) + f(1) \quad (1)$$
$$x = \frac{t}{a}, \quad y = at : \left(\frac{t}{a} + \frac{a}{t}\right) f(at) = f(t^2) + f(a^2) \quad (2)$$
$$x = a^2 t, \quad y = t : \left(a^2 t + \frac{1}{a^2 t}\right) f(t) = f(a^2 t^2) + f\left(\frac{1}{a^2}\right) \quad (3)$$
$$x = y = at : \left(at + \frac{1}{at}\right) f(at) = f(a^2t^2) + f(1) \quad (4)$$
In order to eliminate $f(t^2)$, take the difference of Eq. (1) and Eq. (2); from Eq. (3) and Eq. (4) eliminate $f(a^2t^2)$; then by taking a linear combination, eliminate $f(at)$ as well:
$$\begin{aligned} & \left(t + \frac{1}{t}\right) f(t) - \left(\frac{t}{a} + \frac{a}{t}\right) f(at) = f(1) - f(a^2) \text{ and} \\ & \left(a^2 t + \frac{1}{a^2 t}\right) f(t) - \left(at + \frac{1}{at}\right) f(at) = f\left(\frac{1}{a^2}\right) - f(1), \text{ so} \\ & \left(\left(at + \frac{1}{at}\right)\left(t + \frac{1}{t}\right) - \left(\frac{t}{a} + \frac{a}{t}\right)\left(a^2 t + \frac{1}{a^2 t}\right)\right) f(t) \\ & = \left(at + \frac{1}{at}\right)\left(f(1) - f(a^2)\right) - \left(\frac{t}{a} + \frac{a}{t}\right)\left(f\left(\frac{1}{a^2}\right) - f(1)\right). \end{aligned}$$

Notice that on the left-hand side, the coefficient of $f(t)$ is nonzero and does not depend on $t$:
$$\left( at + \frac{1}{at} \right) \left( t + \frac{1}{t} \right) - \left( \frac{t}{a} + \frac{a}{t} \right) \left( a^2 t + \frac{1}{a^2 t} \right) = a + \frac{1}{a} - \left( a^3 + \frac{1}{a^3} \right) < 0.$$
After dividing by this fixed number, we get:
$$f(t) = C_1 t + \frac{C_2}{t} \qquad (5)$$
where the numbers $C_1$ and $C_2$ are expressed in terms of $a$, $f(1)$, $f(a^2)$ and $f(1/a^2)$, and they do not depend on $t$.
The functions of the form Eq. (5) satisfy the equation:
$$\begin{align*} \left(x + \frac{1}{x}\right) f(y) &= \left(x + \frac{1}{x}\right) \left(C_1 y + \frac{C_2}{y}\right) \\ &= \left(C_1 xy + \frac{C_2}{xy}\right) + \left(C_1 \frac{y}{x} + C_2 \frac{x}{y}\right) \\ &= f(xy) + f\left(\frac{y}{x}\right). \end{align*}$$

Solution 2. We start with an observation. If we substitute $x = a \neq 1$ and $y = a^n$ in the assumption, we obtain
$$f(a^{n+1}) - \left(a + \frac{1}{a}\right)f(a^n) + f(a^{n-1}) = 0.$$
For the sequence $z_n = a^n$, this is a homogenous linear recurrence of the second order, and its characteristic polynomial is
$$t^2 - \left(a + \frac{1}{a}\right)t + 1 = (t - a)(t - \frac{1}{a})$$
with two distinct nonzero roots, namely $a$ and $1/a$. As is well-known, the general solution is
$$z_n = C_1 a^n + C_2 \left(\frac{1}{a}\right)^n$$
where the index $n$ can be as well positive as negative. Of course, the numbers $C_1$ and $C_2$ may depend of the choice of $a$, so in fact we have two functions, $C_1$ and $C_2$, such that
$$f(a^n) = C_1(a) \cdot a^n + \frac{C_2(a)}{a^n} \text{ for every } a \neq 1 \text{ and every integer } n. \quad (1)$$
The relation Eq. (1) can be easily extended to rational values of $n$, so we may conjecture that $C_1$ and $C_2$ are constants, and whence
$$f(t) = C_1 t + \frac{C_2}{t}.$$
As it was seen in the previous solution, such function indeed satisfy the assumption.

a solution of the assumption. In order to make our formulas simpler,
define
$$f_0(c) = f(x) - f(1) \cdot x.$$
This function is another one satisfying the assumption and the extra
constraint $f_0(1) = 0$. Repeating the same argument on linear recur-
rences, we can write
$$f_0(a) = K(a)a^n + \frac{L(a)}{a^n}$$
with some functions $K$ and $L$. By substituting $n = 0$, we can see that
$$K(a) + L(a) = f_0(1) = 0 \text{ for every } a.$$
Hence,
$$f_0(a^n) = K(a)\left(a^n - \frac{1}{a^n}\right).$$
Now take two numbers $a > b > 1$ arbitrarily and substitute
$x = (a/b)^n$ and $y = (ab)^n$ in the assumption:
$$\left(\frac{a^n}{b^n} + \frac{b^n}{a^n}\right) f_0((ab)^n) = f_0(a^{2n}) + f_0(b^{2n}), \text{ so}$$
$$\left(\frac{a^n}{b^n} + \frac{b^n}{a^n}\right) K(ab) \left((ab)^n - \frac{1}{(ab)^n}\right) = K(a) \left(a^{2n} - \frac{1}{a^{2n}}\right) + K(b) \left(b^{2n} - \frac{1}{b^{2n}}\right),$$
or equivalently
$$K(ab)\left(a^{2n} - \frac{1}{a^{2n}} + b^{2n} - \frac{1}{b^{2n}}\right) = K(a)\left(a^{2n} - \frac{1}{a^{2n}}\right) + K(b)\left(b^{2n} - \frac{1}{b^{2n}}\right). \quad (2)$$

$$K(ab) = K(a).$$
Then Eq. (2) reduces to $K(a) = K(b)$. Hence, $K(a) = K(b)$ for all $a > b > 1$.
Fix $a > 1$. For every $x > 0$ there is some $b$ and an integer $n$ such that $1 < b < a$ and $x = b^n$. Then
$$f_0(x) = f_0(b^n) = K(b)\left(b^n - \frac{1}{b^n}\right) = K(a)\left(x - \frac{1}{x}\right).$$
Hence, we have
$$f(x) = f_0(x) + f(1)x = C_1x + \frac{C_2}{x}$$
with
$$C_1 = K(a) + f(1) \text{ and } C_2 = -K(a).$$