5
$\begingroup$

This is a PDE taken from a Maple document. Mathematica DSolve currently unable to solve it.

I wanted to verify Maple solution using NDSolve. This is string of length 1, fixed on the left, and free to move on the right. Given an initial position and let go.

Here is the specs of the PDE

Solve for $0<x<1, t>0$ the wave PDE $$ -u_{tt} + u(x,t)= u_{xx} + 2 e^{-t} \left( x - \frac{1}{2} x^2 + \frac{1}{2} t - 1 \right) $$

With boundary condition

\begin{align*} u(0,t) &= 0 \\ \frac{\partial u(1,t)}{\partial x} &= 0 \end{align*}

And initial conditions

\begin{align*} u(x,0) &= x^2-2 x \\ u(x,1)&= u\left(x,\frac{1}{2}\right) + e^{-1} \left( \frac{1}{2} x^2-x\right) - \left( \frac{3}{4} x^2- \frac{3}{2}x \right) e^{\frac{-1}{2}} \end{align*}

The tricky part in this, is that no initial velocity is given. But only initial position at $t=0$, and then a relation on the solution at 2 different times is give instead.

NDSolve complain with that dreaded error

Boundary condition is not specified on a single edge of the boundary of the computational domain.

And I do not know how to get rid of it. Here is the code

ClearAll[u, x, t];
pde = -D[u[x, t], {t, 2}] + u[x, t] == 
           D[u[x, t], {x, 2}] + 2*Exp[-t]*(x - (1/2)*x^2 + (1/2)*t - 1);

bc = {u[0, t] == 0, Derivative[1, 0][u][1, t] == 0};

ic = {u[x, 0] == x^2 - 2*x, 
    u[x, 1] == u[x, 1/2] + ((1/2)*x^2 - x)*Exp[-1] - ((3*x^2)/4 - (3/2)*x)* Exp[-2^(-1)]};

sol = NDSolve[{pde, ic, bc}, u, {x, 0, 1}, {t, 0, 1}]

Here is the Maple code and the analytical solution it gives

pde := -diff(u(x, t), t, t) + u(x, t) = 
       diff(u(x, t), x, x)+ 2*exp(-t)*(x-(1/2)*x^2+(1/2)*t-1);
ic  := u(x, 0) = x^2-2*x, 
       u(x, 1) = u(x, 1/2)+((1/2)*x^2-x)*exp(-1)-(3/4*(x^2)-3/2*x)*exp(-1/2);
bc  := u(0, t) = 0, eval(diff(u(x, t), x), {x = 1}) = 0;
pdsolve([pde, ic, bc],u(x,t))

$$ u(x,t) = -\frac{e^{-t}}{2} (x^2-2 x) (t-2) $$

Here is animation of Maple solution, which I wanted to verify

mapleSol[x_, t_] := -(Exp[-t]/2) (x^2 - 2 x) (t - 2)
Manipulate[
 Plot[mapleSol[x, t], {x, 0, 1}, PlotRange -> {{0, 1}, {-1, .1}}],
 {{t, 0, "time"}, 0, 10, .1}
 ]

enter image description here

Any suggestion how to get rid of the error from NDSolve?

Using V 12 on windows 10. ps. I solved this by hand also, but can't get Maple solution, and my solution looks wrong. I still need to find out why.

$\endgroup$
  • 1
    $\begingroup$ What Maple delivers here is simply excellent! $\endgroup$ – rmw Jul 23 at 10:08
  • $\begingroup$ @rmw yes. Nice solution. I struggled to find out how to get same solution analytically for few hrs and could not. $\endgroup$ – Nasser Jul 23 at 10:14
5
$\begingroup$

I don't have time to make a solid answer, but this seems to work :

pde = -D[u[x, t], {t, 2}] + u[x, t] == 
   D[u[x, t], {x, 2}] + 2*Exp[-t]*(x - (1/2)*x^2 + (1/2)*t - 1) + 
    NeumannValue[0, x == 1];

bc = {u[0, t] == 0};

ic = {u[x, 0] == x^2 - 2*x
   , PeriodicBoundaryCondition[
    u[x, t] - (((1/2)*x^2 - x)*Exp[-1] - ((3*x^2)/4 - (3/2)*x)*
        Exp[-2^(-1)])
    , t == 1 && 0 < x < 1
    , Function[xy, xy - {0, 1/2}]]};

U = NDSolveValue[{pde, ic, bc}, u, {x, 0, 1}, {t, 0, 1}];
Plot3D[U[x, t], {x, 0, 1}, {t, 0, 1}, AxesLabel -> {x, t, u}]  

enter image description here

Plot[{
  U[x, 0]
  , U[x, 1/2]
  , U[x, 1]
  , U[x, 1/2] + ((1/2)*x^2 - x)*Exp[-1] - ((3*x^2)/4 - (3/2)*x)*
    Exp[-2^(-1)]}, {x, 0, 1}, 
 PlotStyle -> {Red, Green, Directive[Blue, AbsoluteThickness[7]], 
   Directive[Black, Dashed, AbsoluteThickness[3]]}, 
 PlotLegends -> "Expressions"]  

enter image description here

Here is the error :

Plot3D[Evaluate[-D[U[x, t], {t, 2}] + 
   U[x, t] - (D[U[x, t], {x, 2}] + 
     2*Exp[-t]*(x - (1/2)*x^2 + (1/2)*t - 1))], {t, 0, 1}, {x, 0, 1}]  

enter image description here

The method automaticaly chosen by NDSolve is Method -> {"PDEDiscretization" -> {"FiniteElement"}} (as opposed to Method] -> {"PDEDiscretization" -> {"MethodOfLines", "SpatialDiscretization" -> {"FiniteElement", femopts}}}) . This is the reason why one can impose boundaries condition on the variable "time".

Note also that the term "PeriodicBoundaryCondition" is a little bit misleading because the source of the "boundary condition" does not need to be a boundary.

$\endgroup$
  • 1
    $\begingroup$ What a tricky solution! Thanks. That means in a rectangular region it is possible to define coupled boundary conditions which aren't periodic at all. For example PeriodicBoundaryCondition[u[x, t] - f[x,t] , t == 1 && 0 < x < 1 , Function[xy, xy - {x0, t0}]] is somthing like u[x,1]+f[x,1]==u[x-x0,1-t0]+f[ x-x0,t-t0]??? $\endgroup$ – Ulrich Neumann Jul 23 at 8:48
  • $\begingroup$ Nice answer (+1). @UlrichNeumann, you should have a look at the ref page of PeriodicBoundaryCondition That explains a bit better what PBC does and has examples. $\endgroup$ – user21 Jul 23 at 11:15
  • $\begingroup$ @user21 Thank you for your hint, I tried to understand the documentation. Do you think my conclusion (see comment) is ok? $\endgroup$ – Ulrich Neumann Jul 23 at 11:25
  • $\begingroup$ @UlrichNeumann, I think that is correct but I hesitate a bit because I have not tried it. $\endgroup$ – user21 Jul 23 at 12:07
  • $\begingroup$ @user21 Thanks, I also try to build an example... $\endgroup$ – Ulrich Neumann Jul 23 at 12:10
1
$\begingroup$

Just an extended comment:

If you change the second bc to NeumannValue Mathematica is able to solve the modified initial value problem u[x, 1] ==(* u[x,1/2]+*) ((1/2)*x^2 - x)*Exp[-1] - ((3*x^2)/4 - (3/2)*x)*Exp[-2^(-1)]

pde = -D[u[x, t], {t, 2}] + u[x, t] ==D[u[x, t], {x, 2}] + 2*Exp[-t]*(x - (1/2)*x^2 +(1/2)*t - 1) +NeumannValue[0, x == 1];

bc = {u[0, t] == 0};

ic = {u[x, 0] == x^2 - 2*x, 
u[x, 1] ==(* u[x,1/2]+*) ((1/2)*x^2 - x)*Exp[-1] - ((3*x^2)/4 - (3/2)*x)*Exp[-2^(-1)]};

U = NDSolveValue[{pde, ic, bc}, u, {x, 0, 1}, {t, 0, 1} ];
Plot3D[U[x, t], {x, 0, 1}, {t, 0, 1}, AxesLabel -> {x, t, u}]

enter image description here

The coupling u[x,1],u[x,1/2] still remains unsolved!

$\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.