calculus - Conceptual question on substitution in integration

In calculus we learn about the substitution method of integrals, but I haven't been able to prove that it works. I mainly don't see how manipulations of differentials is justified, i.e how $dy/dx = f(x)$ means that $dy = dx * f(x)$ so $dx * f(x)$ can be substituted for $dy$, since I thought that $dy/dx$ is merely notation and that $dy$ and $dx$ don't actually exist.

Answer

Essentially we want to show:

$$\int_a^bf(y)dy=\int_{g^{-1}(a)}^{g^{-1}(b)}f(g(x))g'(x)dx$$
for a strictly increasing continuous function $g$ (if $g$ is decreasing, we take $-g$ and absorb the minus sign into swapping the limits). If $\{y_0,\ldots,y_n\}$ is a partition of $[a,b]$, $y_{j-1}\leq y_j^*\leq y_j$, then a Riemann sum for the left integral is
$$\sum_{j=1}^nf(y^*_j)(y_j-y_{j-1}).$$
Put $y_j:=g(x_j)$. So $\{x_0,\ldots,x_n\}$ is a partition of $[g^{-1}(a),g^{-1}(b)]$. By the mean value theorem, there exists $x_j^*\in[x_{j-1},x_j]$ such that
$$g'(x_j^0)=\frac{g(x_j)-g(x_{j-1})}{x_j-x_{j-1}}\iff y_j-y_{j-1}=g'(x_j^*)(x_j-x_{j-1}).$$
Now remember that all Riemann sums converge to the integral (by definition of a function being Riemann integrable), so we may choose $y_j^*$ so that $x_j^*=g^{-1}(y_j^*)$. Hence we have
$$\sum_{j=1}^nf(y^*_j)(y_j-y_{j-1})=\sum_{j=1}^nf(g(x_j^*)g'(x_j^*)(x_j-x_{j-1}).$$

Now we simply recognise that the right hand side is a Riemann sum for the right integral and that $\max|y_j-y_{j-1}|\to0$ as $\max|x_j-x_{j-1}|\to0$.

Sorry for such a lengthy answer - the basic takeaway is that you can prove it, and the notation we use is just a shorthand. While $dy$ and $dx$ are not actual objects as you rightly point out, they act similarly enough in a lot of ways that we often treat them as such.