calculus - Why can we treat infinitesimals as real numbers in integration by substitution?

During integration by substitution we normally treat infinitesimals as real numbers, though I have been made aware that they are not real numbers but merely symbolic, and yet we still can, apparently, treat them as real numbers. For instance, consider we want to integrate the expression $3x(x^4+1)^3$. A common way to do this is to let $u=x^4+1$, where $\frac{du}{dx}=4x^3$, and thus $du=4x^3dx$ which is appropriately used in our substitution to obtain $\int3x(u)^4 du$, and then we simply directly integrate this new integrand. However, while I understand the process and why we do it in such a manner, I am perplexed as to why we can still rigorously treat the infinitesimals as real numbers. So, my question is if anyone can elaborate on exactly why it is logically rigorous to treat infinitesimals as real numbers during substitution for integration.

(Note: My question does not concern as to what "dx" means in integration simply because my question is defined in the prospect of treating infinitesimal derivatives as ratios specifically in integration by substitution, where other questions do not specifically address. )

Answer

The issue is that there is a whole bunch of measure theory hidden behind those "legitimate" calculations. If you just learned how to integrate elementary functions, you have too much to learn before you get there.

The key in this is the Radon-Nikodym Theorem, which says that if you have two measures $\mu,\nu$ on a measurable space $X$ such that $\mu$ is absolutely continuous with respect to $\nu$, then there exists a unique density function (which we denote by $\frac{d\mu}{d\nu}$) such that for any function $g$ integrable with respect to both $\mu$ and $\nu$,
$$
\int_X g \, d\mu = \int_X g \frac{d\mu}{d\nu} \, d\nu.
$$
The function $\frac{d\mu}{d\nu} : X \to \mathbb R$ is called the Radon-Nikodym derivative of $\mu$ with respect to $\nu$.

So in other words, in the rigorous treatment, the infinitesimals are not treated as infinitely small quantities, but rather as measures; while doing a change of variables, it can be shown that the Radon-Nikodym derivative you obtain is in fact the derivative of the function you use to change variables (i.e. in $du = f'(x) dx$, $\, f'(x) = \frac{du}{dx}$ where $du$ and $dx$ are two measures on the real line and $u = f(x)$ ; $f$ must be a diffeomorphism).

Now let me assume that this was not satisfactory for you; there is also the field of mathematics called non-standard analysis, which defined hyperreal numbers; I suggest you look at the Wikipedia page on this for more details. It allows the treatment of infinitesimals as quantities so you can play around with them.

Hope that helps,