The complete derivation is slightly long and mathematical in nature so I'll just mention some key points.
Before Planck came into the picture, physicists wanted to study light-matter interaction. To do so, they experimentally realised a nearly ideal blackbody (using what experimentalists call a 'Jean Cube') and obtained the graph of spectral radiance vs wavelength. The experimentalists were happy, but the theorists needed a framework to understand the graphs obtained.
Instead of using spectral radiance, they found it more useful to start off with energy density function $\rho (f)$. Now, to find the energy density function, the approach the physicists took was:
- Count the number of waves that can fit inside Jean's cube. Because radiation inside Jean's cube gets reflected back and forth, it was reasonable to assume the waves inside Jean's cube were stationary waves.
- Find the average energy of these waves in thermal equilibrium. This is where statistical mechanics comes to help us. The physicists assumed the entities to be the standing waves here and not particles.
- Using the results from 1) and 2), we can calculate the energy between the frequency range $f$ and $f+df$. This effectively gives us the expression of the energy density function $\rho (f)$.
Nice. Step 1 is calculated through a few arguments but nothing to worry about here. What makes the difference is step 2.
Before Planck, everyone assumed energy to be continuous. Therefore, we calculate the average energy of the system (using statistical mechanics) as:
$$
\bar{E}=\frac{\int_0^{\infty} E e^{-\frac{E}{k T}} d E}{\int_0^{\infty} e^{-\frac{E}{k T}} d E}
$$
It can be easily shown that this becomes $\bar{E} = kT $.
Using the fact: total energy within the interval $df$ equals $($total number of waves in $df$ interval$)$ times $($average energy per wave$)$ times $($volume of jeans cube$)^{-1}$ and moving from a frequency domain to wavelength domain, we end up obtaining:
$$\rho(\lambda)d\lambda = \frac{8 \pi kT} {\lambda ^{4}} d\lambda$$
(don't worry about how we exactly got this right now and all what's needed is just the result we obtain)
One can convert the energy density function to spectral radiance, and then it turns out that our theoretical model shows $\frac{1}{\lambda ^{4}}$ behaviour and this is why we obtain the 'Ultraviolet catastrophe'. Clearly, for shorter wavelengths - around the ultraviolet range - the spectral radiance seems to diverge.
This was a big issue during those times. This meant that either our theoretical models were wrong or incomplete. Max Planck claimed that the theoretical models were incomplete in the sense that for long wavelengths, our theoretical graphs tend to match with experimental graphs (so we aren't completely wrong). Because for long wavelengths we get some overlap, it is fair to assume $\bar{E} \approx kT$ is correct in this range. Max Planck successfully solved this puzzle by assuming energy to be discrete. Nobody knows what thought process made Planck make the following assumption (but it worked somehow):
$$E=n\varepsilon; n \in N$$
(where $N$ is the set of natural numbers)
Because the energy is now discrete, this means, going back to step 2, our calculations become:
$$
\bar{E}=\frac{\int_0^{\infty} E e^{-\frac{E}{k T}} d E}{\int_0^{\infty} e^{-\frac{E}{k T}} d E} \Rightarrow \bar{E}=\frac{\sum_{n=0}^{\infty} n \varepsilon e^{-\frac{n \varepsilon}{k T}}}{\sum_{n=0}^{\infty} e^{-\frac{n \varepsilon}{k T}}}
$$
The latter can be calculated using techniques from differentiation and geometric progression to obtain:
$$
\bar{E}=\frac{\varepsilon}{e^{\frac{\varepsilon}{k T}}-1}
$$
Okay, we got that by assuming energy is discrete. That's fine but what's important are the following observations:
- From experimental graphs, as frequency goes to zero, $\bar{E}$ tends to $kT$ and when frequency tends to infinity, $\bar{E}$ tends to $0$.
- From the above equation $\bar{E}=\frac{\varepsilon}{e^{\frac{\varepsilon}{k T}}-1}$, we find that as $\varepsilon$ goes to zero, $\bar{E}$ tends to $kT$ and when $\varepsilon$ goes to infinity, $\bar{E}$ tends to $0$.
This means that $\varepsilon$ and frequency $f$ are somehow related, so we postulate the following relation:
$$ \varepsilon = hf$$
where $h$ is just some constant.
Okay, we use $\bar{E}=\frac{\varepsilon}{e^{\frac{\varepsilon}{k T}}-1}$ and proceed what we did earlier and we obtain:
$$
\rho(\lambda) d \lambda=\frac{8 \pi h c}{\lambda^5} \frac{d \lambda}{\left[e^{\frac{h c}{\lambda k T}}-1\right]}
$$
The behaviour of the energy density function is now very much similar to the behaviour seen in the experimental data. All that's left to do is calibrate the constant $h$ so that the theoretical graph matches the experimental data exactly. Doing so gave $h$ a value that we today associate with Planck's constant.
Supplementing this along with Nadav's answer to this question effectively should explain everything as to why quantisation solves the ultraviolet catastrophe.