Inhomogeneous dynamic Euler-Bernoulli beam equation with discontinuous parameters

Question

The code

x0 = 0.25;  T = 20;   u1 = -0.03;  u2 = 0.07;  u3 = -0.04;
a = 1/100;   t0 = 5;  omega = 2;
a = 0.01;  dis[x_] := a/(Pi (x^2 + a^2))
P[t_] := If[t <= t0, Sin[omega t], 0]
u[t_] := u1 HeavisideTheta[t - 0.8] + 
  u2 HeavisideTheta[t - 1.64] + u3 HeavisideTheta[t - 3.33]

pde = a D[w[x, t], {x, 4}] + D[w[x, t], {t, 2}] - 
   P[t] dis[x - x0];
sol = NDSolve[{pde == 0, w[0, t] == u[t], w[1, t] == 0, 
    Derivative[2, 0][w][0, t] == 0, Derivative[2, 0][w][1, t] == 0, 
    w[x, 0] == 0, Derivative[0, 1][w][x, 0] == 0}, 
   w[x, t], {x, 0, 1}, {t, 0, 80}, Method -> "StiffnessSwitching"];

gives an error

NDSolve::eerr: Warning: scaled local spatial error estimate of 246.5944594961422` at t = 80.` in the direction of independent variable x is much greater than the prescribed error tolerance. Grid spacing with 25 points may be too large to achieve the desired accuracy or precision. A singularity may have formed or a smaller grid spacing can be specified using the MaxStepSize or MinPoints method options.

I think it's because of the boundary conditions on derivatives. Have tried Automatic, MethodOfLines, etc., does not help. Tried

Method -> {"MethodOfLines", 
      "SpatialDiscretization" -> {"TensorProductGrid", 
        "MinPoints" -> 100}}

works fine with second order equation subjected to bc containing only first order derivative. Any thoughts, hints?

Answer 1

eerr is a warning, not an error, it just suggests the possibility of trouble and doesn't always mean the output you obtained is wrong. Indeed, the solution given by NDSolve with the default setting seems to be erroneous, but according to my test, with a spatial grid dense enough e.g. 50 (BTW it seems to be better to use a even number), NDSolve won't give too bad a solution, even if the warning is still there:

appro = With[{k = 1000}, ArcTan[k #]/Pi + 1/2 &];
unitStepExpand = Simplify`PWToUnitStep@PiecewiseExpand@# &;

x0 = 25/100; T = 20; u1 = -3/100; u2 = 7/100; u3 = -4/100;
a = 1/100; t0 = 5; omega = 2;
a = 1/100; dis[x_] := a/(Pi (x^2 + a^2))
P[t_] = unitStepExpand@If[t <= t0, Sin[omega t], 0];
u[t_] = u1 HeavisideTheta[t - 8/10] + u2 HeavisideTheta[t - 164/100] + 
    u3 HeavisideTheta[t - 333/100] /. HeavisideTheta -> UnitStep;

pde = a D[w[x, t], {x, 4}] + D[w[x, t], {t, 2}] - P[t] dis[x - x0];

(* I decrease the value of tend 
   because the solution seems to come to steady state long before 80. *)
tend = 10;

mol[n_Integer, o_:"Pseudospectral"] := {"MethodOfLines", 
  "SpatialDiscretization" -> {"TensorProductGrid", "MaxPoints" -> n, 
    "MinPoints" -> n, "DifferenceOrder" -> o}}

fsol = NDSolveValue[{pde == 0, w[0, t] == u[t], w[1, t] == 0, 
      Derivative[2, 0][w][0, t] == 0, Derivative[2, 0][w][1, t] == 0, w[x, 0] == 0, 
      Derivative[0, 1][w][x, 0] == 0} /. UnitStep -> appro, w, {x, 0, 1}, {t, 0, tend}, 
    MaxSteps -> Infinity, Method -> mol[50, 4]]; // AbsoluteTiming

Plot[fsol[0, t], {t, 0, tend}, PlotRange -> .05]
Plot3D[fsol[x, t], {x, 0, 1}, {t, 0, tend}, PlotRange -> 2, PlotPoints -> 40, 
 ColorFunction -> "AvocadoColors", Lighting -> {{"Ambient", White}}]
Manipulate[Plot[fsol[x, t], {x, 0, 1}, PlotRange -> 2], {t, 0, tend}]

I modfied the parameters with the appro and unitStepExapnd because of the reason mentioned here.

Well, I admit there still exists slight yet suspicious oscillation in the solution, and obtaining a highly accurate numerical solution for the problem seems to be pending, so let's wait and see if someone will come up with a better approach.

Answer 2

I think that maybe Mathematica can't solve the boundary condition on the second order derivative problem.(I'm not sure, waiting for other expert's answer!) So, I tried to rewrite your code by defining s[x, t] == D[w[x, t], {x, 2}]

Therefore, your pde and bc becomes:

pde = {a D[s[x, t],{x, 2}]+D[w[x, t],{t, 2}]-P[t] dis[x- x0]==0,
       s[x, t] == D[w[x, t], {x, 2}]};
bc = {w[0, t] == u[t], w[1, t] == 0, s[0, t] == 0, s[1, t] == 0, 
      w[x, 0] == 0, Derivative[0, 1][w][x, 0] == 0};
sol = NDSolve[{pde, bc}, {s, w}, {x, 0, 1}, {t, 0, 80}];
Plot3D[w[x, t] /. sol, {x, 0, 1}, {t, 0, 80}]

with result: