# Inhomogeneous dynamic Euler-Bernoulli beam equation with discontinuous parameters

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?

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 = SimplifyPWToUnitStep@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.

• Thank you very much! Yes, you are right, the oscillations are very suspicious. – Aso Oct 29 '16 at 16:51

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: • I have also tried this, but it does not help. – Aso Oct 29 '16 at 12:56
• "I think that maybe Mathematica can't solve the boundary condition on the second order derivative problem." Well, I don't think so, here's a problem that is quite similar to OP's. – xzczd Oct 29 '16 at 14:31
• BTW, which version do you use? I just tested the code on Cloud but can't reproduce your result, only got the (0,1) CoefficientArrays::poly:` warning. – xzczd Oct 29 '16 at 16:15
• I use Wolfram Mathematics 10.1. – Aso Oct 29 '16 at 16:51
• Thanks for your helpful information, xzczd! I use version 10.0.1.0. – tablecircle Oct 30 '16 at 1:26