The strange thing about the maximum in Planck's law I read that it makes a difference whether you calculate $\frac{dE(\lambda) }{d \lambda}=0$ or $\frac{dE(\omega)}{d \omega}=0$ in the sense that the maximum energy density with respect to the wavelength does not coincide with the frequency maximum that one would assume to be at $\omega_{max} = \frac{c}{\lambda_{max}}$. Actually, I do not understand why this is so ( Now, I am only interested in a pure physical explanation, mathematically this is clear and due to the chain rule). Does anybody know how to explain this odd thing?
Somehow I feel that the core of the question has not been completely answered. Although I can easily look up wikipedia, where the need for the chain rule is explained, I am rather interested in understanding where the following argument breaks down, which seems to be why this question has caused some confusion in the past(as you can see by googling this question):
So we have $E(\omega)$ the energy radiated at a given frequency. Now this function has a maximum somewhere, so there is a frequency where a maximal amount of energy is emitted. In other words: If you add the sum of the photon's energies at each frequency that are emitted, you will notice that the maximum is reached at this frquency. Now $E(\lambda)$ tells you basically the same for the wavelength, but again: We know where at which frequency the maximal amount of energy is radiated, so we know the corresponding wavelength. 
I think there is a need to explain this. 
If anything is unclear, please let me know.
 A: Let's revisit what the term "spectral energy density" means.It means the amount of energy emitted in a infinitesimal $d \lambda$ or $d \nu$.
Now due the the relationship  $\lambda=\frac{c}{\nu}$,we can find that $d \lambda=-\frac{c}{\nu^2}d \nu$.
The $\nu^2$ in the denominator leads to the given phenomenon.We can think of problem as finding the box($d \lambda$) which has the maximum height in the Planck curve.
Obviously the orange box is the required maximum.However if we were to represent the given graph in terms of the frequency,the size of boxes($d \lambda$ or $d \nu$) would change due to the $\nu^2$ in the denominator.As a result the box which corresponded to the maximum height in the wavelength curve may not correspond to the max. height in the frequency curve.This is exactly what happens and the graphs peak at two different points.
A: There is no much Physics involved. The two graphs below ($\sin(x)$ and $\sin(x^{-1})/x^2$) summarize it. The largest areas are the green and the blue. In the first graph, the areas have equal widths. In the second graph, the areas are reversed in order and the widths are not the same.

A: Here's an answer using no math other than counting, which is all we need to explain this phenomenon.
Images generated with this script.
Throw a pair of dice 1000 times. The individual results look like this:
10,6,11,4,6,8,7,5,7,5,11,6,9,8,11,5,7,11,7,6,8,9,8,6,9,8,9,7,6,12...

Now if we plot the frequency of each result, we observe that seven turns up most often:

The reason is that the result seven can be achieved in many different ways (1+6, 2+5, 3+4, 4+3, 5+2 and 6+1) whereas e.g. two can only be obtained when both dice give exactly 1.
One might now expect that plotting the reciprocals of these results (ranging from $\tfrac1{12}$ to $\tfrac12$) would show a peak at $\tfrac17 \approx 0.143$. Well, turns out this does not happen – the peak is actually rather at $0.12$:

The reason for this is binning: because the reciprocals fall more densely around low values (the spacing between $\tfrac1{12}$ and $\tfrac1{11}$ is only $0.0076$, while the spacing between $\tfrac13$ and $\tfrac12$ is $0.16$). Hence we actually get more different results in the $0.12$ bin than in the $0.14$ bin.
Now, for this dice example we could remove this issue by making so many so small bins that every possible result has a bin of its own:

In this view, the peak is at $0.14$. Most of the bins are zero.
But this only works because the individual results are hard-quantised: you can't get two dice to sum to e.g. $8.734$.
A physical distribution like the Planck spectrum is continuous, i.e. no matter how small we make the bins, we will always have “different events” in each bin. And setting up bins with equal distribution in frequency space gives a different result from equal-spaced bins in wavelength space.
A: The physics is simple: the wavelengh is not proportional to the frequency, but a "distorted" function of the latter, so the density functions (factors at $d\nu$ and $d\lambda$) are different.
A: Yes it confuses many people. In that sense, we should have 2 wien's law; one giving the maximum in wavelength and the other in frequency, plus a Wien's law which will be defined with some new variable such as $\nu^2$ or $\lambda^3$. Everytime, we do that, we will obtain a new wien's law. Physics which depend so much on our choice of variables should be just abandoned because it produces confusion. There was a nice publication about it, where they propose a new definition
https://arxiv.org/pdf/1109.3822.pdf
A: I'll try to point by point address the question.
So we have $E(\omega)$ the energy radiated at a given frequency.
Actually, no, we have an energy per unit frequency.  It's like a density.  For a mass density $\rho$ you have a mass per unit volume, and the density can vary from place to place, so to get the total mass you break up space into a bunch of little volumes, each of which has density constant to within some tolerance, then you take that density and that volume and multiply them together to get the mass of that piece, and add up all the pieces, or $m=\iiint\rho dx dy dz$.
Similarly you have energy $U$ and energy density $u$.  But they come in multiple kinds, you can have an energy per unit wavelength $u_\lambda$ or an energy per unit frequency $u_\nu$, and these two different energy densities (per wavelength and per frequency) don't even have the same units.  For an energy per wavelength $u_\lambda$, you take all the possible wavelengths, chop the interval into pieces, each of which the density $u_\lambda$ doesn't change much and multiply the density $u_\lambda$ by the interval of wavelengths in the pieces $\Delta \lambda$ to get $u_\lambda \Delta \lambda$ the total energy $dU$ in that piece, then you add that up for every piece of the total interval to get the total energy:
$$U=\int_0^\infty u_\lambda d\lambda.$$
Similarly for frequency, you take the total frequency range $(0,\infty)$, chop the interval into pieces, each of which the density $u_\nu$ doesn't change much and multiply the density $u_\nu$ by the interval of frequencies in the piece $\Delta \nu$ to get $u_\nu \Delta \nu$ the total energy in that piece, then you add that up for every piece $dU$ of the total interval to get the total energy:
$$U=\int_0^\infty u_\nu d\nu.$$
You could even do this in reverse.  Take the interval of all frequencies and break it into pieces, then for each piece take the total energy in that piece $dU$, and divide it by the size of the piece $\Delta \nu$ to get $dU/\Delta\nu$, and notice that when $\Delta \nu$ is small enough, the quotient $dU/\Delta \nu$ doesn't change anymore ($dU$ gets smaller too when you consider a smaller range of frequencies).  That quotient is the energy density $u_\nu$.
If you do this, you see that for any piece of the total interval, there is a real amount of energy in the interval, $dU$, and that there is a real lower frequency(high wavelength) and a real higher frequency (low wavelength) and everything you think is fine.  The only difference is the length of the interval of wavelengths is different than the length of the interval of frequencies.
For instance if you went from 1Hz to 2Hz in frequency, then $\Delta \nu=$2Hz-1Hz=1Hz, so you get $u_\nu = dU/\Delta \nu=dU/$1Hz.  For the exact same interval of frequencies there is an interval of wavelengths that covers the same piece of the total.  Since $\lambda \nu=c$ (your equation was off by a factor of $2\pi$ since you used angular frequency), we have a lower wavelength of $c/2Hz$ and an upper wavelength of $c/1Hz$, these are wavelengths and have units of wavelength, and in this case it is about $1.5\times10^8 m$.  The same $dU$ is divided by 1Hz to get $u_\nu$ but is divided by $1.5\times10^8 m$ to get $u_\lambda$, doesn't seem like a big deal.  However, what if your frequency range was from $2 Hz$ to $3 Hz$ , then there is a new $dU$ for this range.  And the $\Delta \nu$ is again $3Hz-2Hz=1Hz$, but the $\Delta \lambda$ equals $c/2Hz-c/3Hz$ which is about $5\times10^7 m$.  So the exact same bit of energy $dU$ is divided by a smaller range of wavelengths, so for the exact same interval of physical waves $u_\lambda$ is bigger than $u_\nu$ (bigger than it was for a slightly different interval that in frequency space was just as long).  The energy involved $dU$ is exactly the same always, and the cut off wavelengths obey $\lambda \nu = c$, but the $\Delta \lambda$ and the $\Delta \nu$ are different from interval to interval.  So the densities are different, since the density is the ratio of $dU$ the real energy to $\Delta \lambda$ (or $\Delta \nu$) for energy per wavelength (or energy per frequency respectively).
Now this function has a maximum somewhere, so there is a frequency where a maximal amount of energy is emitted. 
The function is not energy at a frequency, it is the energy $dU$ divided by the length $\Delta \nu$ of the frequency interval in question.  It isn't the energy.  It is a maximum when $dU$ gets larger than any other $dU$ for intervals with the same sized $\Delta \nu$.  Imagine taking a very very small interval $\Delta \nu$, and moving it to different places (but keeping the same size) until you find a place where $dU$ for that interval is largest.  That's what it means for a density to be at its maximum.  That's the frequency where a window of frequencies of fixed frequency spread gives the most energy.  And these windows of fixed frequency spread (say 1Hz, or one mHz or one $\mu$Hz) don't have a fixed wavelength spread as you move them to different places, so it has nothing to do with where density per wavelength is maximum.
To find where $u_\lambda$ is maximal, note that is the energy $dU$ divided by the length $\Delta \lambda$ of the frequency interval in question.  It isn't the energy.  It is a maximum when $dU$ gets larger than any other $dU$ for intervals with the same sized $\Delta \lambda$.  Imagine taking a very very small interval $\Delta \lambda$, and moving it to different places (but keeping the same size) until you find a place where $dU$ for that interval is largest.  That's what it means for a density to be at its maximum.  That's the wavelength where a window of wavelengths of fixed wavelength spread gives the most energy.  And these windows of fixed wavelength spread (say 1mm, or one m or one km) don't have a fixed frequency spread as you move them to different places, so it has nothing to do with where density per frequency is maximum.
In other words: If you add the sum of the photon's energies at each frequency that are emitted, you will notice that the maximum is reached at this [frequency].
No, and this time it is wrong on many additional levels, this brings up photons, and energy density depends on total energy in the interval, it is not a per photon measure, that would be a totally different issue.  And we aren't adding up energy per frequency, we are adding up the energy from intervals of frequencies, the result for each interval is $dU=u_\nu\Delta\nu$.
Now $E(\lambda)$ tells you basically the same for the wavelength
Except if we use the same intervals, they now have the same $dU$ for each interval, but the $\Delta \lambda$ for each interval is different, so $u_\lambda$ is different.
We know where at which frequency the maximal amount of energy is radiated, so we know the corresponding wavelength.
We might know which interval $dU$ gave the most energy, and if all the intervals were of fixed size $\Delta \nu$ in frequency space then indeed we can take $u_\nu=dU/\Delta\nu$ and the largest one will be where $u_\nu$ is largest.  But to get $u_\lambda$, we have to compute $u_\lambda=dU/\Delta\lambda$ and the $\Delta \lambda$ are different for these intervals, so the one that made $u_\nu=dU/\Delta\nu$ largest is not the one that makes $u_\lambda=dU/\Delta\lambda$ largest.
It really comes down to knowing what an energy density is, an energy divided by an amount of stuff (spread of wavelengths, spread of frequency, region of volume, etc.)
The bounty was for well documented answers, so I tried to include all the details so that this answer by itself answers everything.  The only thing you really need to know is what a density is, which is something you put into an integal, energy per wavelength $u_\lambda$ is something that gives
$$U=\int u_\lambda d\lambda,$$
for the total energy for any range of wavelengths.  And energy per frequency $u_\nu$ is something that gives
$$U=\int u_\nu d\nu,$$
for the total energy for any range of frequencies.  Absolutely everything else follows from knowing that, which is just the definition of a density.
A: There are many different choices for the variable $x$ describing the kind of plane harmonic EM wave (its frequency), and each may lead to different function $I_x(x)$ that is to express the original spectral distribution of the studied radiation. If new variable $x' = T(x)$ is introduced, this means that generally the equality
$$
I_x(x) = I_{x'}(x')~~~(1)
$$
need not be satisfied. If it was, the maximum of $I_x$ would have $x$ corresponding to $x'$ that maximizes $I_{x'}$; but it often is not.
This is because the transformation to other spectral functions is done rather by requiring that
$$
I_x(x) dx= I_{x'}(x')dx',
$$
together with using prescription for the value of
$$
\frac{dx'}{dx},
$$
which follows from the relation $T(x)$ defining $x'$ based on $x$.
Since these two functions $I_x, I_{x'}$ do not need to satisfy (1), their maxima may correspond to different kinds of EM wave, and then obviously the wave whose wavelength maximizes $I_{\lambda}(\lambda)$ need not be the same wave whose frequency maximizes $I_f(f)$.
