This post contains several code blocks, you can copy them easily with the help of importCode.
Analytic Solution
The analytic solution can be obtained with LaplaceTransform and FourierSinCoefficient. First, make a Laplace transform on the equation and b.c.s and plug in the i.c.s:
Clear[f];
f[x_] = x (1 - x);
eqn = D[u[t, x], {t, 2}] + D[u[t, x], {x, 4}] == 0;
ic = {u[0, x] == f@x, Derivative[1, 0][u][0, x] == 0};
bc = {u[t, 0] == 0, u[t, 1] == 0, Derivative[0, 2][u][t, 0] == 0,
Derivative[0, 2][u][t, 1] == 0};
teqn = LaplaceTransform[{eqn, bc}, t, s] /. Rule @@@ ic
Now we have an ODE, solve it with DSolve:
tsol = u[t, x] /. First@DSolve[teqn/.
HoldPattern@LaplaceTransform[a_, __] :> a, u[t, x], x] // Simplify
Notice the replacement HoldPattern@LaplaceTransform[a_, __] :> a is necessary because DSolve has trouble in handling expression containing LaplaceTransform. The last step is to transform the solution back, but sadly InverseLaplaceTransform can't handle tsol. At this point, one work-around is to turn to numeric inverse Laplace transform, you can use this or this package for this task. But for your specific problem, we can circumvent the issue by expanding tsol with Fourier sine series:
easyFourierSinCoefficient[expr_, t_, {a_, b_}, n_] :=
FourierSinCoefficient[expr /. t -> t + a, t, n,
FourierParameters -> {1, Pi/(b - a)}] /. t -> t - a
easyTerm[t_, {a_, b_}, n_] := Sin[Pi/(b - a) n (t - a)]
term = easyTerm[x, {0, 1}, n];
coe = easyFourierSinCoefficient[tsol, x, {0, 1}, n]
$$-\left(i\left(\frac{(1+i) (-1)^n e^{i \sqrt{2} \sqrt{s}}}{(1+i) \pi n+i \sqrt{2} \sqrt{s}}\right.\right....$$
coe still looks complex, but inspired by those (-1)^ns in it, we split it to odd and even part and simplify:
oddcoe =
Simplify[coe /. n -> 2 m - 1, m > 0 && m ∈ Integers] /. m -> (1 + n)/2
(* (8 s)/(n^3 π^3 (n^4 π^4 + s^2)) *)
evencoe = Simplify[coe /. n -> 2 m, m ∈ Integers] /. m -> n/2
(* 0 *)
InverseLaplaceTransform can handle the series form of the transformed solution without difficulty:
soloddterm = Function[n, #] &@InverseLaplaceTransform[oddcoe term, s, t]
(* Function[n, (8 Cos[n^2 π^2 t] Sin[n π x])/(n^3 π^3)] *)
To find the final solution, just summate:
solgenerator[n_] := Compile[{t, x}, #] &@Total@soloddterm@Range[1, n, 2];
sol = solgenerator[200];
Animate[Plot[sol[t, x], {x, 0, 1}, PlotRange -> .3], {t, 0, 1}]
The animation is similar to the one in the subsequent solution so I'd like to omit it.
Fully NDSolve-based Numeric Solution
Go back to the old-fashioned "TensorProductGrid", set "DifferentiateBoundaryConditions" -> {True, "ScaleFactor" -> 100} (or NDSolve will set "ScaleFactor" to 0 so the inconsistent b.c.s will be severely ignored, for more information, check the obscure tutorial) and DifferenceOrder -> 2:
mol[n_Integer, o_:"Pseudospectral"] := {"MethodOfLines",
"SpatialDiscretization" -> {"TensorProductGrid", "MaxPoints" -> n,
"MinPoints" -> n, "DifferenceOrder" -> o}}
mol[tf:False|True,sf_:Automatic]:={"MethodOfLines",
"DifferentiateBoundaryConditions"->{tf,"ScaleFactor"->sf}}
tu = NDSolveValue[{eqn, ic, bc}, u, {t, 0, 10}, {x, 0, 1},
Method -> Union[mol[27, 2], mol[True, 100], {Method -> StiffnessSwitching}], MaxSteps -> Infinity];
Animate[Plot[tu[t, x], {x, 0, 1}, PlotRange -> .3], {t, 0, 10}]
NDSolve spits out the NDSolveValue::eerr warning, but in many cases NDSolveValue::eerr isn't a big deal, and the result indeed looks OK:
Remark
The {Method -> StiffnessSwitching} option is not necessary in v9, but becomes necessary to obtain a reasonable solution at least since
v12.3.
The performance of this approach backslides severely at least since v12.3. For {t, 0, 2}, the timing is about 16 seconds in v9, but about 275 seconds in
v12.3. I haven't found out the reason so far, but the posteriori error estimation is 604.34 in v12.3, which is smaller than the 1802.02
in v9.0.1.
Partly NDSolve-based Numeric Solution
Theoretically we can also set "DifferentiateBoundaryConditions" -> False to avoid the inconsistent b.c.s being ignored, but strangely NDSolve spits out the icfail warning and fails. I'm not sure about reason, but found that we can manually discretize the spatial derivative and solve the obtained ODE set with NDSolve to avoid the issue.
First, let's define a function pdetoode that discretizes PDEs to ODEs (Additionally, though not related to OP's problem, I've also define a function pdetoae that discretizes differential equations to algebraic equations based on pdetoode. A diffbc function is defined to transform b.c. to (almost) equivalent ODE, in certain cases you may need it to avoid the fragile DAE solver of NDSolve. A rebuild function is also created to combine the list of InterpolatingFunctions to a single InterpolatingFunction):
Clear[fdd, pdetoode, tooderule, pdetoae, diffbc, rebuild]
fdd[{}, grid_, value_, order_, periodic_] := value;
fdd[a__] := NDSolve`FiniteDifferenceDerivative@a;
pdetoode[funcvalue_List, rest__] :=
pdetoode[(Alternatives @@ Head /@ funcvalue) @@ funcvalue[[1]], rest];
pdetoode[{func__}[var__], rest__] := pdetoode[Alternatives[func][var], rest];
pdetoode[front__, grid_?VectorQ, o_Integer, periodic_: False] :=
pdetoode[front, {grid}, o, periodic];
pdetoode[func_[var__], time_, {grid : {__} ..}, o_Integer,
periodic : True | False | {(True | False) ..} : False] :=
With[{pos = Position[{var}, time][[1, 1]]},
With[{bound = #[[{1, -1}]] & /@ {grid}, pat = Repeated[_, {pos - 1}],
spacevar = Alternatives @@ Delete[{var}, pos]},
With[{coordtoindex =
Function[coord,
MapThread[Piecewise[{{1, PossibleZeroQ[# - #2[[1]]]},
{-1, PossibleZeroQ[# - #2[[-1]]]}}, All] &, {coord, bound}]]},
tooderule@Flatten@{
((u : func) | Derivative[dx1 : pat, dt_, dx2___][(u : func)])[x1 : pat, t_,
x2___] :> (Sow@coordtoindex@{x1, x2};
fdd[{dx1, dx2}, {grid}, Outer[Derivative[dt][u@##]@t &, grid],
"DifferenceOrder" -> o, PeriodicInterpolation -> periodic]),
inde : spacevar :>
With[{i = Position[spacevar, inde][[1, 1]]}, Outer[Slot@i &, grid]]}]]];
tooderule[rule_][pde_List] := tooderule[rule] /@ pde;
tooderule[rule_]@Equal[a_, b_] :=
Equal[tooderule[rule][a - b], 0] //. eqn : HoldPattern@Equal[_, _] :> Thread@eqn;
tooderule[rule_][expr_] := #[[Sequence @@ #2[[1, 1]]]] & @@ Reap[expr /. rule]
pdetoae[funcvalue_List, rest__] :=
pdetoae[(Alternatives @@ Head /@ funcvalue) @@ funcvalue[[1]], rest];
pdetoae[{func__}[var__], rest__] := pdetoae[Alternatives[func][var], rest];
pdetoae[func_[var__], rest__] :=
Module[{t},
Function[pde, #[
pde /. {Derivative[d__][u : func][inde__] :>
Derivative[d, 0][u][inde, t], (u : func)[inde__] :> u[inde, t]}] /. (u : func)[
i__][t] :> u[i]] &@pdetoode[func[var, t], t, rest]]
diffbc[rst__][a : _List | _Equal] := diffbc[rst] /@ a
diffbc[dvar : {t_, order_} | (t_) .., sf_: 0][a_] /; sf =!= t := sf a + D[a, dvar]
rebuild[funcarray_, grid_?VectorQ, timeposition_: 1] :=
rebuild[funcarray, {grid}, timeposition]
rebuild[funcarray_, grid_, timeposition_?Negative] :=
rebuild[funcarray, grid, Range[Length@grid + 1][[timeposition]]]
rebuild[funcarray_, grid_, timeposition_: 1] /; Dimensions@funcarray === Length /@ grid :=
With[{depth = Length@grid},
ListInterpolation[
Transpose[Map[Developer`ToPackedArray@#["ValuesOnGrid"] &, #, {depth}],
Append[Delete[Range[depth + 1], timeposition], timeposition]],
Insert[grid, Flatten[#][[1]]["Coordinates"][[1]], timeposition]] &@funcarray]
The syntax of pdetoode is as follows: 1st argument is the function to be discretized (which can be a list i.e. pdetoode can handle PDE system), 2nd argument is the independent variable in the resulting ODE system (usually it's the variable playing the role of "time" in the underlying model), 3rd argument is the list of spatial grid, 4th argument is difference order, 5th argument is to determine whether periodic b.c. should be set or not. (5th argument is optional, the default setting is False. )
Notice pdetoode is a general purpose function. You may feel some part of the source code confusing. To understand it, just notice the following truth:
a /. a | b[m_] :> {m} outputs {}.Derivative[][u] outputs u.
Then discretize eqn, ic and bc and remove redundant equations:
lb = 0; rb = 1;
(* Difference order of x: *)
xdifforder = 2;
points = 25;
grid = Array[# &, points, {lb, rb}];
(* There're 4 b.c.s, so we need to remove 4 equations from every PDE/i.c.,
usually the difference equations that are the "closest" ones to the b.c.s
are to be removed: *)
removeredundant = #[[3 ;; -3]] &;
(* Use pdetoode to generate a "function" that discretizes the spatial derivatives of
PDE(s) and corresponding i.c.(s) and b.c.(s): *)
ptoofunc = pdetoode[u[t, x], t, grid, xdifforder];
odeqn = eqn // ptoofunc // removeredundant;
odeic = removeredundant/@ptoofunc@ic;
odebc = bc // ptoofunc;
sollst = NDSolveValue[{odebc, odeic, odeqn}, u /@ grid, {t, 0, 10}, MaxSteps -> Infinity];
(* Rebuild the solution for the PDE from the solution for the ODE set: *)
sol = rebuild[sollst, grid];
Animate[Plot[sol[t, x], {x, 0, 1}, PlotRange -> .3], {t, 0, 10}]
The animation is similar to the one in the aforementioned solution so I'd like to omit it. This approach seems to be more robust than the fully NDSolve-based one, because even if the xordereqn i.e. the difference order for spatial derivative is set to 4, it's still stable, while the fully NDSolve-based one becomes wild when t is large.
