What is the probability that taking at random $a$ and $b$ $\leq m$, $a$ divides $b$? While doing some computational experiments with recurrences, I came with this simple question: if we take random natural numbers $a$ and $b$ $\leq m \in \mathbb{N}$, what is probability that the recurrence $ax_{n}-bx_{n-1}=0, x_{0}=1$ is an integral sequence? (i.e $x_{n} \in \mathbb{Z}\;\forall\;n)$. Since the recurrence is trivial ($x_{n}=\left(\frac{b}{a}\right)^{n}$), the sequence is integral if $a|b$.
So the fundamental question is,

What is the probability that taking at random natural numbers $a$ and $b$ $\leq m$, $a$ divides $b$?

It is know that the probability that a random number $a$, less or equal than $b$, is coprime with $b$ is $\frac{\phi(b)}{b}$, where $\phi(b)$ is the Euler's totient function but I don't know how to use this argument for the previous question.
 A: Presumably your natural numbers exclude $0$.
If the number of divisors of $n$ is $d(n)$ then the probability is $$\frac1{m^2}\sum\limits_{n=1}^{m} d(n)$$
which is close to $$\frac1m\left(\log_e(m)+2\gamma-1\right)$$ where $\gamma\approx 0.5772156649$ is the Euler–Mascheroni constant. This has a limit of $0$ as $m$ increases without bound.
How close?  This is essentially Dirichlet's divisor problem divided by $m^2$ so the difference is at worst $O\left(m^{-1.685}\right)$  in Big-O notation
A: Your sample space takes $(a,b)$ from $\{1,\ldots,m\}^2$. Assuming you are selecting uniformly, the $m^2$ outcomes are equally likely.
Now count how many outcomes have $a\mid b$. If $a=1$, there are $m$ values for $b$. If $a=2$, there are $\lfloor m/2\rfloor$ values for $b$. And so on.
So the probability is $$\frac{1}{m^2}\sum_{j=1}^{m}\left\lfloor\frac{m}{j}\right\rfloor$$
A: $\newcommand{\floor}[1]{\lfloor #1 \rfloor}$
We can choose $a$ and $b$ from the set $(1, 2, 3, \dots, m)$. Thus, there are $m^2$ ways to choose them.

If $b$ has to be divisible by $a$, we can do casework based on what $a$ is. When $a=1$, $b$ can be anything from $1$ to $m$, giving $m$ values possible. When $a=2$, then there are $\floor{\frac{m}{2}}$ ways to choose $b$. Similarly, when $a=n$, there are $\floor{\frac{m}{n}}$ ways to choose $b$ out of the set $(1, 2, 3, \dots, m)$.

So, our final answer is $\frac{\floor{\frac{m}{2}} + \floor{\frac{m}{3}} + \floor{\frac{m}{4}} + \dots + \floor{\frac{m}{m}}}{m^2}$, which is in summation form $$\boxed{\frac{1}{m^2} \sum_{n=1}^{m} \floor{\frac{m}{n}}}$$ $\blacksquare$
A: A numerical exercise to get an intuition with Monte Carlo simulations:
import numpy as np
import matplotlib.pyplot as plt
from tqdm import tqdm

def iteration(m:int):
    arlenght = 1000*m
    draws = np.random.randint(low=1, high=m+1, size=(2*arlenght,1))
    A = draws[:arlenght]
    B = draws[arlenght:]
    return(np.count_nonzero(A%B == 0)/arlenght)
m_max = 1000
probs = []
for m in tqdm(range(2,1000)):
    probs.append(iteration(m))

We get:

Following @Henry's answer, the comparison with the approximation using the Euler-Mascheroni constant:
def analytical_approx(n):
    return (np.log(n) + 2*np.euler_gamma - 1)/n
approx = np.array([analytical_approx(m) for m in range(2,m_max+1)])

The comparison yields:

And the ratios seem to converge in mean to one:

