Honestly, the argument you're making here is a mess - the question is based on bad premises. So let me show you how to do it properly, and hopefully that will resolve your confusion.
You're right that in order for a wavefunction to be normalized, it must satisfy
$$\int_\text{all space} P(x)\mathrm{d}x
= \int_\text{all space} \psi^*(x)\psi(x)\mathrm{d}x = 1\tag{1}$$
But this statement:
On top of that, any $\psi$ can also be expressed as $\psi\psi^∗$
is not correct. Given a function $\psi(x)$, you can write $\psi(x)\psi^*(x)$, but that's a different function.
Anyway, given that your wavefunction can be written
$$\psi(x) = a\phi_1(x) + b\phi_2(x)$$
then you just plug that into the normalization condition (1) and get
$$\int_0^L \bigl(a^* \phi_1^*(x) + b^* \phi_2^*(x)\bigr)\bigl(a \phi_1(x) + b \phi_2(x)\bigr)\mathrm{d}x = 1$$
which expands to
$$\begin{multline}
a^*a \int_0^L \phi_1^*(x)\phi_1(x)\mathrm{d}x
+ a^*b \int_0^L \phi_1^*(x)\phi_2(x)\mathrm{d}x \\
+ b^*a \int_0^L \phi_2^*(x)\phi_1(x)\mathrm{d}x
+ b^*b \int_0^L \phi_2^*(x)\phi_2(x)\mathrm{d}x
= 1
\end{multline}\tag{2}$$
Now you can use the identity
$$\int_0^L\phi_1^*(x)\phi_2(x)\mathrm{d}x = \int_0^L\phi_2^*(x)\phi_1(x)\mathrm{d}x = 0$$
which follows from the fact that $\phi_1$ and $\phi_2$ are orthogonal functions (it's not enough that they are eigenfunctions of an operator, they have to be orthogonal), and the identity
$$\int_0^L\phi_1^*(x)\phi_1(x)\mathrm{d}x = \int_0^L\phi_2^*(x)\phi_2(x)\mathrm{d}x = 1$$
which simply reflects the fact that $\phi_1$ and $\phi_2$ are normalized. (Check for yourself that this is the same as the normalization condition, equation (1).) With these two identities, equation (2) reduces to
$$\lvert a\rvert^2 + \lvert b\rvert^2 = 1$$
The conceptual question I had was that if we have the probability squared here, is it that or the square root of that probability that is your normalization constant?
That all depends, how do you define your normalization constant? It depends on what you're normalizing and how exactly you express it as a function. However you do it, the end requirement for normalization is just that $\lvert a\rvert^2 + \lvert b\rvert^2 = 1$.
As far as using the specific sinusoidal form for the $\phi_i$, you can do that in this case, because you're given enough information to figure out that the eigenfunctions are in fact sinusoidal. But you don't really need to know that they are sinusoidal for the preceding argument to work; all you need to know is that the $\phi_i$s are orthonormal.
So here is the abstract approach:
$$ \langle \psi | H | \psi \rangle = \frac{1}{5}\bigg( \langle \phi_1 | H | \phi_1 \rangle + 2\langle \phi_1 | H | \phi_2 \rangle + 2\langle \phi_2 | H | \phi_1 \rangle + 4 \langle \phi_2 | H | \phi_2 \rangle \bigg) \,.$$
Now you know that $H|\phi_1\rangle = E_1 |\phi_1 \rangle$ and $H|\phi_2\rangle = E_2 |\phi_2 \rangle$ --- or rather, you can easily check that the functions you've given are indeed eigenstates of the Hamiltonian:
$$ \frac{-\hbar^2}{2m} \frac{d^2}{dx^2} \phi_n = \frac{n^2 \pi^2 \hbar^2}{2m L^2} \phi_n \equiv E_n \phi_n \,.$$
The functions $\phi_n$ are also normalised, as you can check, and are orthogonal to one another --- this must be the case, because they are the eigenfunctions of a Hermitian operator (with different eigenvalues). Hence the expression above becomes:
$$\langle \psi | H | \psi \rangle = \frac{1}{5}\bigg( E_1\langle \phi_1 | \phi_1 \rangle + 2E_2 \langle \phi_1 | \phi_2 \rangle + 2E_2\langle \phi_2 | \phi_1 \rangle + 4E_2 \langle \phi_2 | \phi_2 \rangle \bigg) \,.$$
$$\langle \psi | H | \psi \rangle = \frac{1}{5}\bigg( E_1 + 4E_2\bigg) \,.$$
Substituting in the form of the energies gives:
$$\langle \psi | H | \psi \rangle = \frac{17}{5} \frac{\pi^2 \hbar^2}{2m L^2} \,.$$
This is, to me, easier than computing the integral you've given, although the integral you've given is correct (or almost --- the Hamiltonian should have factor of $\hbar$ squared in front of the second derivative). If you tried to compute the integral, you would find a good deal of cancellation due to the orthogonality of the functions involved. If you're good at quickly spotting when an integral vanishes, e.g.:
$$ \int_0^L \sqrt{\frac{2}{L}} \sin \left(\frac{\pi x}{L }\right) \sqrt{\frac{2}{L}} \sin \left(\frac{2\pi x}{L }\right) \,dx = 0$$
then the above prescription might not seem any simpler. But as far as cleanness of approach goes, it's nicer to invoke the orthogonality of the eigenfunctions --- which is a result of central importance in this kind of problem, and one which you will prove in any introductory QM course --- before delving into explicit computations.
================================================================================
NB: I anticipate that you may not have met the above notation yet. There isn't anything to it. For our purposes, just take $\langle \psi | H |\psi \rangle$ to mean
$$ \int \psi^*(x) H \psi(x) \,dx \,,$$
from which you should be able to see how the first line follows. The statement that two functions are orthogonal amounts to $\langle \phi_1 | \phi_2 \rangle = 0$, whilst the statement that a function is normalised amounts to $\langle \phi_1 | \phi_1 \rangle = 1$.
Best Answer
The integrals you have to compute are not that difficult. First of all, the spatial parts of your wavefunctions are real (this is wise and always possible in one dimension), and therefore complex conjugates just do not matter. So, the remaining integral you have on the rhs is (up to some expression dependent on t (S(t)) which we will leave for now):
$$ I(t) = S(t) \int_0^a\Phi_n(x)\Phi_m(x) dx = S(t) \frac{2}{a} \int_0^a \operatorname{sin}\Big( \frac{n\pi x}{a} \Big) \operatorname{sin}\Big( \frac{m\pi x}{a} \Big)dx$$
Now, you have to make a substitution $ t = \frac{\pi x}{a} $ to get
$$ I(t) = S(t) \frac{2}{\pi} \int_0^\pi \operatorname{sin} ( n x ) \operatorname{sin} ( m x )dx$$
First of all, if either n or m are 0 then I(t) = 0 and your initial state is normalized.
Now you integrate by parts (or by mathematica if you can't by parts) to obtain
$$ I(t) = S(t) \frac{2}{\pi} [ -(\frac{1}{m} \sin{n x} \cos{m x}) |_0^\pi + \frac{n}{m} \int_0^\pi \cos{n x} \cos({m x}) dx]$$
The first term is 0 (cause sin is 0 in all multiplicities of $\pi$). Now we can proceed with integration by parts of the second term.
$$ I(t) = S(t) \frac{2}{\pi} ( \frac{n}{m} \int_0^\pi \cos{n x} \cos({m x}) dx) = S(t) \frac{2}{\pi} [ (\frac{n}{m^2} \cos{n x} \sin{m x}) |_0^\pi $$ $$ + \frac{n^2}{m^2} \int_0^\pi \sin{n x} \sin({m x}) dx] = \frac{n^2}{m^2} I(t)$$
So, there are two possibilities - either n=m or I(t)=0 (n=-m is excluded as negative labels are not in the basis because sin(-nx) = -sin(nx) - linear dependency). What you get is ORTHOGONALITY RELATION - if eigenvectors have different labels (e.g. 1 and 2 as in your case) they are orthogonal in $L^2([0,a])$ which happens to be your Hilbert space. There are theorems which state that this is always true for Hamiltonian eigenstates - corresponding theory of differental operators is known as Strum-Liouville theory (for finite dimensional Hilbert spaces this is a trivial property of self-adjoined operators).
Now, to the second part of your question. The modulus squared of wavefunction is by definition probability density of finding a particle in an interval $dx$
So, the probability of finding particle in [0,a/2] is simply
$$\int_0^{a/2} |f(x,t)|^2 dx = \frac{1}{2}\int_0^{a/2} |\Phi_1(x)|^2 + |\Phi_2(x)|^2 dx + $$
$$ \frac{1}{2}\int_0^{a/2} \Phi_1(x) \Phi_2(x) \operatorname{exp} \Big( (-1^2+2^2)\frac{-i\pi^2 h t}{2ma^2} \Big)+ \Phi_2(x) \Phi_1(x) \operatorname{exp} \Big( (-2^2+1^2) \frac{-i\pi^2 h t}{2ma^2} \Big) $$ $$ = \frac{1}{2}\int_0^{a/2} |\Phi_1(x)|^2 + |\Phi_2(x)|^2 dx + \cos{\frac{3\pi^2 h t}{2ma}} \int_0^{a/2} \Phi_1(x) \Phi_2(x) dx $$
First integral is simply 1 because we take half of the area of $sin(x)^2$ and $sin(2x)$, and add it. Second integral can be easily computed using double 'by parts' formula from above and replacing $\pi$ with $\pi/2$ in integral limits. One gets $\frac{4}{3 \pi}$ Therefore finally
$$\int_0^{a/2} |f(x,t)|^2 dx = \frac{1}{2}+\frac{4}{3 \pi} \cos{\frac{3\pi^2 h t}{2ma}} $$
As $\frac{4}{3 \pi} < \frac{1}{2}$ the result makes sense. As the body is in mixture of two states, the probability is no longer constant in time.