1
$\begingroup$

I'd like to solve a partial differential equation on a infinite strip. My code is the following:

eqn = D[u[x, t], {t, 2}] + 2 D[u[x, t], {t, 1}] + u[x, t] == 
   D[u[x, t], {x, 2}];
bc = {Derivative[1, 0][u][0, t] == 0, 
   Derivative[1, 0][u][2 π, t] == 0};
ic = {u[x, 0] == Cos[x], Derivative[0, 1][u][x, 0] == Abs[Cos[x/2]]};
Ω = Rectangle[{0, 2 π}, {0, Infinity}];
DSolve[{eqn, bc, ic}, u[x, t], {x, t} ∈ Ω]

Mathematica returns no result -- what could have possible gone wrong?

$\endgroup$
  • $\begingroup$ Rectangle[{0, 2 \[Pi]}, {0, Infinity}] makes no sense for Mathematica. In particular, Graphics[Rectangle[{0, 2 \[Pi]}, {0, Infinity}]] produces an error message. $\endgroup$ – user64494 May 8 at 8:17
  • 1
    $\begingroup$ Related: mathematica.stackexchange.com/q/189706/1871 $\endgroup$ – xzczd May 9 at 5:57
2
$\begingroup$

Infinity does not work with Rectangle. A proper definition can just use Reals instead of Infinity with Elements. Rectangle is not needed, used by Wolfram Inc. for that purpose.

DSolve has on itself an options section for domains.

As shown in the Basic Example for DSolve it is preferable to solve for a pure function u instead of a function u[x,t].

A simple check of the boundary and initial conditions shows, they are inconsistent. The initial condition are inconsistent too. Derivative[0,1] is the derivative to time. It is set to a time-independent function in x with frequency one half.

eqn = D[u[x, t], {t, 2}] + 2 D[u[x, t], {t, 1}] + u[x, t] == 
   D[u[x, t], {x, 2}];
bc = {Derivative[1, 0][u][0, t] == 0, 
   Derivative[1, 0][u][2 π, t] == 0};
ic = {u[x, 0] == Cos[x], Derivative[0, 1][u][x, 0] == Abs[Cos[x/2]]};
ndsol = NDSolve[{eqn, bc, ic}, u, {x, 0, 2 π}, {t, 0, 4.5}]

Mathematica does solve this with the message NDSolve::ibcinc. error message

Plot3D[Evaluate[u[x, t] /. First@%63], {x, 0, 2 Pi}, {t, 0, 4.5}, 
 PlotRange -> All, AxesLabel -> Automatic]

Graphics representation

Boundary conditions are matched well:

bc[[All, 1]] /. t -> 0 /. First@ndsol

(* {-0.000171153, 0.000171153} *)

The boundary values are constant the value for all time.

Initial conditions too:

ic[[All, 1]] /. x -> 0 /. First@ndsol
{1., 1.}

The solution prefers one over the other initial condition on the complete domain.

For example if all bc and is conditions are set to zero, then the message disappears.

Everything is OK with the ic and bc for the solution and it looks:

Graphics

If it is not possible to get on other paths to the correct bc and ic. Fit them to the general solution of the partial differential equation.

| improve this answer | |
$\endgroup$
4
$\begingroup$

DSolve (version 8.0) can't find a solution. Do it numerically.

ndsol = NDSolve[{eqn, bc, ic}, u, {x, 0, 2 Pi}, {t, 0, 10}, 
   MaxStepSize -> {5*10^-3, 10^-2}, AccuracyGoal -> 5, 
   PrecisionGoal -> 5]

Plot3D[Evaluate[u[x, t] /. First@ndsol], {x, 0, 2 Pi}, {t, 0, 10}, 
   PlotRange -> All]

enter image description here

Test

bc[[All, 1]] /. t -> 0 /. First@ndsol

(*   {-5.1341*10^-13, 5.23173*10^-13}   *)

ic[[All, 1]] /. x -> 0 /. First@ndsol

(*   {1., 1.}   *)

Plot3D[Evaluate@Chop[eqn[[1]] - eqn[[2]] /. First@ndsol], 
{x, 0, 2 Pi}, {t, 0, 10}, PlotRange -> 10^-5]
| improve this answer | |
$\endgroup$
  • $\begingroup$ Actually you don't need to make the grid that dense, MaxStepSize -> {0.2, Automatic} is enough. $\endgroup$ – xzczd May 9 at 5:21
3
$\begingroup$

Though DSolve can't handle the problem at the moment, the symbolic solution can be found with the help of finite Fourier cosine transform and its inversion:

eqn = D[u[x, t], {t, 2}] + 2 D[u[x, t], {t, 1}] + u[x, t] == D[u[x, t], {x, 2}];
bc = {Derivative[1, 0][u][0, t] == 0, Derivative[1, 0][u][2 π, t] == 0};
ic = {u[x, 0] == Cos[x], 
   Derivative[0, 1][u][x, 0] == 
    Simplify`PWToUnitStep@PiecewiseExpand[Abs[Cos[x/2]], Reals]};

(* Definition of finiteFourierCosTransform isn't included in this post,
   please find it in the link above. *)
tsetlst = Assuming[{n > 0 && n != 2}, 
      Simplify@finiteFourierCosTransform[{eqn, ic}, {x, 0, 2 Pi}, #]] /. 
     Rule @@@ bc & /@ {0, 2, n} /. HoldPattern@finiteFourierCosTransform[a_, __] :> a

tsollst = u[x, t] /. First@DSolve[#, u[x, t], t] & /@ tsetlst

tsol = Piecewise[{tsollst // Most, {n == 0, n == 2}} // Transpose, tsollst // Last] // 
  FullSimplify

mid = inverseFiniteFourierCosTransform[tsol, n, {x, 0, 2 Pi}]

sol = (mid[[2]] /. {n, C} -> {n, 2, 2} // ReleaseHold) + mid /. {{n, C} -> {n, 3, C}, 
   HoldPattern@Piecewise[_, a_] :> a}

$$u(x,t) = \frac{2 e^{-t} t}{\pi }+\frac{e^{-t} \cos (x) (4 \sin (t)+3 \pi (\sin (t)+\cos (t)))}{3 \pi }+\sum _{n=3}^\infty \frac{\cos \left(\frac{n x}{2}\right) \left(8 \cos \left(\frac{n \pi }{2}\right) \sin \left(\frac{n t}{2}\right) (\sinh (t)-\cosh (t) )\right)}{n \left(n^2-1 \right) \pi}$$

Let's compare it with the numeric solution:

testfunc = Function[{x, t}, #] &[sol /. C -> 30 // ReleaseHold];

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

nsol = NDSolveValue[{eqn, ic, bc}, u, {t, 0, 4}, {x, 0, 2 Pi}, 
   Method -> mol[50, 4]];

Manipulate[Plot[{testfunc[x, t], nsol[x, t]}, {x, 0, 2 Pi}, PlotRange -> 1, 
  PlotStyle -> {Automatic, {Red, Thick, Dashed}}], {t, 0, 4}]

enter image description here

Remark

  1. Abs is transformed to UnitStep so DSolve can solve faster.

  2. Finite Fourier cosine transform for $n=0$ and $n=2$ cases have been calculated separately because currently finiteFourierCosTransform, which is built on Integrate, can't handle these two special cases properly when they're hidden behind the symbolic n.

  3. sol is just a simplified version of mid.

| improve this answer | |
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.