Include boundary conditions in wave equation solver

Updated answer

I had now little time to look at this again to see why DSolve gave series solution with 1/0 in it. It turned put DSolve made a bug. The sum is correct, but it should be over odd $n$ only and not over all $n$. (Will send bug report to Wolfram support so that they can fix this in next version)

I also wrote the solution by hand to verify that this is the case. So now DSolve gives an exact solution for all $n$ when the sum is manually fixed to be over odd numbers also.

Comparing Fixed DSolve solution with NDSolve

You can obtain direct solution using DSolve.

The trick is the boundary conditions. DSolve has hard time with it starting at not $x=0$. Therefore do a mapping by solving for $z=x+\pi$. This transfers the BC from $-\pi \dots \pi$ to $0 \dots 2\pi$ and now DSolve can indeed solve the PDE exactly.

At the end just replace $z$ back by $x+\pi$.

ClearAll[u, x, t, n, z];
pde = D[u[z, t], {t, 2}] == 4*D[u[z, t], {z, 2}];
ic = {u[z, 0] == 0, Derivative[0, 1][u][z, 0] == (Sin[z])^2};
bc = {u[0, t] == 0, u[2 Pi, t] == 0}; (*notice the shift*)
sol = First@(u[z, t] /. DSolve[{pde, ic, bc}, u[z, t], {z, t}]);
sol = sol /. {K[1] -> n, z -> x + Pi}(*reset the spatial shift*)

The above solution is correct, but sum should be over odd $n$ only. Fixed by hand

MMfixedSolution[x_, t_, max_] := 
 Sum[(16 (-1 + (-1)^n) Sin[n t] Sin[1/2 n(Pi + x)])/(n^2 (-16 + n^2)Pi),{n,1, max, 2}]

The above is the final exact solution. Let us now compare the above solution (using few terms in the Fourier series) with the numerical solution posted above by Henrik Schumacher and plot them one top of each.

ClearAll[u, x, t, n];
weqn = D[u[x, t], {t, 2}] == 4*D[u[x, t], {x, 2}];
ic = {u[x, 0] == 0, Derivative[0, 1][u][x, 0] == Sin[x]^2};
bc = {u[-Pi, t] == 0, u[Pi, t] == 0};
solN = NDSolveValue[{weqn, ic, bc}, u, {x, -Pi, Pi}, {t, 0, 4}, 
   MaxStepSize -> 0.01, 
   Method -> {"PDEDiscretization" -> {"MethodOfLines", 
       "SpatialDiscretization" -> {"TensorProductGrid", 
         "MinPoints" -> 100}}}];

Manipulate[
 Plot[{solN[x, t], MMfixedSolution[x, t, 5]}, {x, -Pi, Pi}, 
  PlotRange -> {{-Pi, Pi}, {-2, 2}}, PlotStyle -> {Red, Blue}, 
  PlotLegends -> {"Numerical", "DSolve solution fixed"}],
 {{t, 0, "time"}, 0, 4, .1, Appearance -> "Labeled"}
 ]

Can't even spot the difference now between the numerical and the exact, and this is only using 1,3,5 (i.e. 3 terms only) in the sum !

Verification of DSolve solution by hand

Solving for $t>0,-\pi<x<\pi$ $$ u_{tt}=c^{2}u_{xx} $$ With BC \begin{align*} u\left( -\pi,t\right) & =0\\ u\left( \pi,t\right) & =0 \end{align*} And initial conditions \begin{align*} u\left( x,0\right) & =0\\ u_{t}\left( x,0\right) & =\sin^{2}\left( x\right) \end{align*}

Let $\xi=x+\pi$. When $x=-\pi,\xi=0$ and when $x=\pi,\xi=2\pi$. In terms of $\xi$, the new pde in $U\left( \xi,t\right) $ becomes

$$ U_{tt}=c^{2}U_{\xi\xi} $$ With BC \begin{align*} U\left( 0,t\right) & =0\\ U\left( 2\pi,t\right) & =0 \end{align*} And initial conditions \begin{align*} U\left( \xi,0\right) & =0\\ U_{t}\left( \xi,0\right) & =\sin^{2}\left( \xi\right) \end{align*}

Let $U=X\left( \xi\right) T\left( t\right) $. The PDE\ becomes \begin{align*} \frac{T^{\prime\prime}X}{c^{2}} & =X^{\prime\prime}T\\ \frac{1}{c^{2}}\frac{T^{\prime\prime}}{T} & =\frac{X^{\prime\prime}} {X}=-\lambda \end{align*} Where $\lambda$ is separation constant. Hence the eigenvalue ODE\ is \begin{align*} X^{\prime\prime}+\lambda X & =0\\ X\left( 0\right) & =0\\ X\left( 2\pi\right) & =0 \end{align*} From the boundary conditions, we see that $\lambda>0$ is the only possible value. Therefore the solution to the above ODE\ is $$ X\left( x\right) =A\cos\left( \sqrt{\lambda}\xi\right) +B\sin\left( \sqrt{\lambda}\xi\right) $$ Since $X\left( 0\right) =0$ then $A=0$ and the solution becomes $X\left( \xi\right) =B\sin\left( \sqrt{\lambda}\xi\right) $. Since $X\left( 2\pi\right) =0$ then for non trivial solution we want $\sqrt{\lambda} 2\pi=n\pi$ or $$ \lambda_{n}=\left( \frac{n}{2}\right) ^{2}\qquad n=1,2,3,\cdots $$ Hence the eigenfunctions are $$ X_{n}\left( \xi\right) =\sin\left( \frac{n}{2}\xi\right) \qquad n=1,2,3,\cdots $$ The time ODE now becomes \begin{align*} T_{n}^{\prime\prime}+c^{2}\lambda_{n}T_{n} & =0\\ T_{n}^{\prime\prime}+c^{2}\left( \frac{n}{2}\right) ^{2}T_{n} & =0\\ T_{n}^{\prime\prime}+\frac{c^{2}n^{2}}{4}T_{n} & =0 \end{align*} Which has the solution $$ T_{n}\left( t\right) =D_{n}\cos\left( \frac{cn}{2}t\right) +E_{n} \sin\left( \frac{cn}{2}t\right) $$ Therefore the complete solution becomes \begin{equation} U\left( \xi,t\right) =\sum_{n=1}^{\infty}\left( D_{n}\cos\left( \frac {cn}{2}t\right) +E_{n}\sin\left( \frac{cn}{2}t\right) \right) \sin\left( \frac{n}{2}\xi\right) \tag{1} \end{equation}

Switching back to $x$ the above becomes

\begin{equation} u\left( x,t\right) =\sum_{n=1}^{\infty}\left( D_{n}\cos\left( \frac{cn} {2}t\right) +E_{n}\sin\left( \frac{cn}{2}t\right) \right) \sin\left( \frac{n}{2}\left( x+\pi\right) \right) \tag{1A} \end{equation}

At $t=0$ the above becomes \begin{align*} 0 & =\sum_{n=1}^{\infty}D_{n}\sin\left( \frac{n}{2}\left( x+\pi\right) \right) \\ D_{n} & =0 \end{align*}

The solution (1A) simplifies to

\begin{equation} u\left( x,t\right) =\sum_{n=1}^{\infty}E_{n}\sin\left( \frac{cn} {2}t\right) \sin\left( \frac{n}{2}\left( x+\pi\right) \right) \tag{2} \end{equation}

Taking time derivative of (2) gives $$ u_{t}\left( x,t\right) =\sum_{n=1}^{\infty}\left( E_{n}\frac{cn}{2} \cos\left( \frac{cn}{2}t\right) \right) \sin\left( \frac{n}{2}\left( x+\pi\right) \right) $$ At $t=0$ the above becomes $$ \sin^{2}\left( x\right) =\sum_{n=1}^{\infty}E_{n}\frac{cn}{2}\sin\left( \frac{n}{2}\left( x+\pi\right) \right) $$ Applying orthogonality gives \begin{align*} \int_{-\pi}^{\pi}\sin^{2}\left( x\right) \sin\left( \frac{n}{2}\left( x+\pi\right) \right) dx & =E_{n}\frac{cn}{2}\int_{-\pi}^{\pi}\sin ^{2}\left( \frac{n}{2}\left( x+\pi\right) \right) dx\\ & =E_{n}\frac{\pi cn}{2} \end{align*}

For the LHS, for $n$ even $\int_{-\pi}^{\pi}\sin^{2}\left( x\right) \sin\left( \frac{n}{2}\left( x+\pi\right) \right) dx=0$. Hence $E_{n}=0$ for all $n$ even. For $n$ odd \begin{align*} \int_{-\pi}^{\pi}\sin^{2}\left( x\right) \sin\left( \frac{n}{2}\left( x+\pi\right) \right) dx & =\frac{16\left( \cos\left( n\pi\right) -1\right) }{\left( n^{2}-16\right) n}\\ & =\frac{16\left( \left( -1\right) ^{n}-1\right) }{\left( n^{2} -16\right) n} \end{align*}

But $n$ is odd, hence the above simplifies more to

$$ \int_{-\pi}^{\pi}\sin^{2}\left( x\right) \sin\left( \frac{n}{2}\left( x+\pi\right) \right) dx=\frac{-32}{\left( n^{2}-16\right) n} $$

Therefore \begin{align} E_{n} & =\frac{2}{n\pi c}\frac{-32}{\left( n^{2}-16\right) n}\nonumber\\ & =\frac{-64}{n^{2}\pi c\left( n^{2}-16\right) }\nonumber \end{align}

Therefore the final solution (2) now becomes $$ u\left( x,t\right) =\sum_{n=1,3,5,\cdots}^{\infty}\frac{-64}{n^{2}\pi c\left( n^{2}-16\right) }\sin\left( \frac{cn}{2}t\right) \sin\left( \frac{n}{2}\left( x+\pi\right) \right) $$

Since $c=2$ then

$$ u\left( x,t\right) =\sum_{n=1,3,5,\cdots}^{\infty}\frac{-32}{n^{2}\pi \left( n^{2}-16\right) }\sin\left( nt\right) \sin\left( \frac{n} {2}\left( x+\pi\right) \right) $$

Which is same as DSolve but sum should be over odd $n$.

You might want to try NDSolveValue instead with the following syntax:

weqn = D[u[x, t], {t, 2}] == 4*D[u[x, t], {x, 2}];
ic = {u[x, 0] == 0, Derivative[0, 1][u][x, 0] == (Sin[x])^2};
bc = {u[-Pi, t] == 0, u[Pi, t] == 0};
sol = NDSolveValue[{weqn, ic, bc}, u, {x, -Pi, Pi}, {t, 0, 4},
   MaxStepSize -> 0.01,
   Method -> {"PDEDiscretization" -> {"MethodOfLines", 
       "SpatialDiscretization" -> {"TensorProductGrid", 
         "MinPoints" -> 100}}}
   ];
Plot3D[sol[x, t], {x, -Pi, Pi}, {t, 0, 4}, Mesh -> None, 
 AxesLabel -> {"x", "t"}]

Assuming symmetric solution u[x,t]==u[-x,t] (particulary Derivative[1,0][u][0,t]) you might, alternatively to @HenrikSchumacher 's answer, use "FiniteElement"

sol = NDSolveValue[{weqn, ic, bc}, u, {x, 0, Pi}, {t, 0, 4}, Method ->{"PDEDiscretization" -> {"MethodOfLines","SpatialDiscretization" -> {"FiniteElement"}}} ];
Plot3D[sol[x, t], {x, 0, Pi}, {t, 0, 4}, Mesh -> None,AxesLabel -> {"x", "t"}]