# Solving a damped wave equation

I am trying to solve the equation

$$\frac{d^2u}{dt^2}-\frac{d^2}{dx^2}\left(c_s^2u+\nu\frac{du}{dt}\right)=0$$

with initial conditions

$$u(x, 0)=0$$

$$\frac{du}{dt}|_{t=0}=0$$

and boundary conditions

$${\frac{du}{dx}}|_{x=0,1}=a\sin{(\omega_d t)}-b\cos{(\omega_d t)}$$

My attempts so far are

ClearAll[u, x, t, a, b, c, w, n];

c = 1;
n = 1;
a = 1;
b = 1;
w = 1;

pde = D[u[t, x], t, t] - c*D[u[t, x], x, x] - n*D[u[t, x], x, x, t] == 0;
ics = {u[0,x]==0};
bcs =
{(D[u[t,x],x] /. x->0) == a*Sin[w*t]-b*Cos[w*t],
(D[u[t,x],x] /. x->1) == a*Sin[w*t]-b*Cos[w*t]};

sol = NDSolveValue[{pde, ics, bcs}, u, {x, 0, 1},{t, 0, 10}]


but I am receiving multiple error messages. Among them:

NDSolveValue::fembdnl: The dependent variable in (u^(0,1))[t,0]==-Cos[t]+Sin[t] in the boundary condition DirichletCondition[(u^(0,1))[t,0]==-Cos[t]+Sin[t],x==0.] needs to be linear.

NDSolveValue::fembdnl: The dependent variable in (u^(0,1))[t,0]==-Cos[t]+Sin[t] in the boundary condition DirichletCondition[(u^(0,1))[t,0]==-Cos[t]+Sin[t],x==0.] needs to be linear.

NDSolveValue::femcmsd: The spatial derivative order of the PDE may not exceed two.

How can I solve it?

• 1. According to your code, The $\rho$s in $\LaTeX$ formula seem to be typos, if so, please correct them, if not, please clarify what's the relationship between $u$ and $\rho$. 2. To find a particular solution, one more initial condition is needed. – xzczd Oct 29 '19 at 12:59
• @xzczd yes,they are typos. Sorry for that, I will correct them. – Alessandro Zunino Oct 30 '19 at 9:48
• With the new added i.c. NDSolve solves the problem without difficulty, ibcinc warning is generated, but it doesn't seem to be a big deal in this case. – xzczd Oct 30 '19 at 10:32
• @xzczd can you provide an example of your working code? Which warning do you get exactly? – Alessandro Zunino Oct 30 '19 at 10:34
• ClearAll[u, x, t, a, b, c, w, n]; c = 1; n = 1; a = 1; b = 1; w = 1; pde = D[u[t, x], t, t] - c*D[u[t, x], x, x] - n*D[u[t, x], x, x, t] == 0; ics = {u[0, x] == 0, D[u[t, x], t] == 0 /. t -> 0}; bcs = {(D[u[t, x], x] /. x -> 0) == a*Sin[w*t] - b*Cos[w*t], (D[u[t, x], x] /. x -> 1) == a*Sin[w*t] - b*Cos[w*t]}; sol = NDSolveValue[{pde, ics, bcs}, u, {x, 0, 1}, {t, 0, 10}]; Plot3D[sol[t, x], {t, 0, 10}, {x, 0, 1}, PlotRange -> All] – xzczd Oct 30 '19 at 10:37

It is necessary to divide the equation into two (second order), add the initial data and boundary conditions. Here is the code without adding boundary conditions

ClearAll[u, x, t, a, b, c, w, n];

c = 1;
n = 1;
a = 1;
b = 1;
w = 1;

pde = {D[u[t, x], t, t] - c*D[u[t, x], x, x] - n*D[v[t, x], t] == 0,
v[t, x] == D[u[t, x], x, x]};
ics = {u[0, x] == 0, v[0, x] == 0, Derivative[1, 0][u][0, x] == 0};
bcs = {(D[u[t, x], x] /. x -> 0) ==
If[t <= 10^-6, 0,
a*Sin[w*t] - b*Cos[w*t]], (D[u[t, x], x] /. x -> 1) ==
If[t <= 10^-6, 0, a*Sin[w*t] - b*Cos[w*t]]};

{U, V} = NDSolveValue[{pde, ics, bcs}, {u, v}, {x, 0, 1}, {t, 0, 10},
Method -> {"IndexReduction" -> Automatic,
"EquationSimplification" -> "Residual",
"PDEDiscretization" -> {"MethodOfLines",
"SpatialDiscretization" -> {"TensorProductGrid",
"MinPoints" -> 141, "MaxPoints" -> 141,
"DifferenceOrder" -> 2}}}]

DensityPlot[U[t, x], {x, 0, 1}, {t, 0, 10},
ColorFunction -> "Rainbow", PlotLegends -> Automatic,
FrameLabel -> Automatic] We can use boundary conditions twice, then there are no warnings

ClearAll[u, x, t, a, b, c, w, n];

c = 1;
n = 1;
a = 1;
b = 1;
w = 1;
f[t_] := If[t < 10^-6, 0, 1];
pde = {D[v[t, x], t] - c*D[u[t, x], x, x] - n*D[v[t, x], x, x] == 0,
v[t, x] == D[u[t, x], t]};
ics = {u[0, x] == 0, v[0, x] == 0};
bcs = {(D[u[t, x], x] /. x -> 0) == (a*Sin[w*t] - b*Cos[w*t]) f[
t], (D[u[t, x], x] /. x -> 1) == a*Sin[w*t] - b*Cos[w*t] f[t]};
bcs1 = {(D[v[t, x], x] /. x -> 0) ==
w (a*Cos[w*t] + b*Sin[w*t]) f[t], (D[v[t, x], x] /. x -> 1) ==
w (a*Cos[w*t] + b*Sin[w*t]) f[t]};

{U, V} = NDSolveValue[{pde, ics, bcs, bcs1}, {u, v}, {x, 0, 1}, {t, 0,
10}, Method -> {"IndexReduction" -> Automatic,
"EquationSimplification" -> "Residual",
"PDEDiscretization" -> {"MethodOfLines",
"SpatialDiscretization" -> {"TensorProductGrid",
"MinPoints" -> 141, "MaxPoints" -> 141,
"DifferenceOrder" -> 2}}}];

{DensityPlot[U[t, x], {x, 0, 1}, {t, 0, 10},
ColorFunction -> "Rainbow", PlotLegends -> Automatic,
FrameLabel -> Automatic, PlotLabel -> "u"],
DensityPlot[V[t, x], {x, 0, 1}, {t, 0, 10},
ColorFunction -> "Rainbow", PlotLegends -> Automatic,
FrameLabel -> Automatic, PlotLabel -> "v"]} • Thank you for your answe, it is very useful. Though, I have a few doubts on your procedure: why did you add an IF condition in the boundary conditions? Is it needed to add the new initial condition v[0, x] == 0 ? Is there a particular reason why you specified all those options to solve the pde? – Alessandro Zunino Oct 27 '19 at 10:25
• The initial conditions Derivative[1, 0][u][0, x] == 0 and boundary conditions are not matched for t = 0. To harmonize these conditions, we excite the boundary conditions immediately at t>0. – Alex Trounev Oct 27 '19 at 12:10
• If I try to run your code with mathematica 12, I receive the error NDSolveValue::bcart: Warning: an insufficient number of boundary conditions have been specified for the direction of independent variable x. Artificial boundary effects may be present in the solution.. Furthermore, using discontinuous boundary conditions apparently causes the solution to have some artifacts at the time point of the discontinuity. – Alessandro Zunino Oct 29 '19 at 9:05
• @AlessandroZunino You can remove If[], but determine Derivative[1, 0][u][0, x]. Otherwise, we have a warning NDSolveValue::ibcinc: Warning: boundary and initial conditions are inconsistent. – Alex Trounev Oct 29 '19 at 9:28
• If I remove the if[] indeed the error you quote appears, even if Derivative[1, 0][u][0, x]==0 is present. Even if this happens, the solutions is still exactly the same you show. Some other errors appear: NDSolveValue::bcart: Warning: an insufficient number of boundary conditions have been specified for the direction of independent variable x. Artificial boundary effects may be present in the solution. – Alessandro Zunino Oct 29 '19 at 10:25