# Strange Boundary NDSolve Problem

Mma10.1 doesn't work with the following boundary problem:

NDSolve[{
y[s, t] == Derivative[2, 0][x][s, t]/(1 + Derivative[1, 0][x][s, t])^(7/2),
x[s, 0] == 0,
Derivative[0, 1][x][s, 0] == 0,
x[0, t] == 0,
y[s, t] == Derivative[0, 2][x][s, t],
y[1, t] == -1 + 5/(1 + Derivative[1, 0][x][1, t])^(5/2)
},
{x[s, t], y[s, t]},
{s, 0, 1}, {t, 0, 1}]


Some messages are printed:

CoefficientArrays::poly: "-(x$11640/(1+x$11639)^(7/2))+y is not a polynomial."


and

NDSolve::femper: "PDE parsing error of {-(x$11640/(1+x$11639)^(7/2))+y,-x$11641+y}. Inconsistent equation dimensions."  • @ WateSoyan : Please clarify your somewhat uncommon way of writing your equations. You should e.g. eliminate y[s,t]. This gives you a non linear version of the wave equation and a strange boundary condition for the second derivative D[x[s,t],{t,2}][1,t]. Ok. If you then consider the linearized wave equation (assuming Abs[D[x[s,t],s]]<<1) you get only the trivial solution x = 0. – Dr. Wolfgang Hintze Jun 2 '15 at 10:42 • I added y as a new varibles since the original equations meet another message: NDSolve::bdord: "Boundary condition 1 + Derivative[0, 2][x][1, t] - 5/(1 + Derivative[1, 0][x][1, t])^(5/2) should have derivatives of order lower than the differential order of the partial differential equation" – WateSoyan Jun 2 '15 at 12:31 ## 2 Answers EDIT Sorry but this development contains an error (replacement of bound2 by bound2a). Thanks to bbgodfrey for notifying me. It turns out that the original bound2 is ill defined in the OP as it contains a derivative of the same order as the differential equation. Only if one makes deliberately the same error as I did (writing dx/dt instaed of dx/ds on the right hand side of bound2), declaring it a modification, a reasonable non trivial result appears. Original text In the first step we write down the mathematical formulation of the problem. Then we transform the equations to Mathematica, and solve the problem. Mathematical formulation The function to be determined is$x(s,t)$. It is governed by this nonlinear partial differential equation $$\text{eq}=\frac{\partial ^2x(s,t)}{\partial t^2}\text{==}\frac{\partial ^2x(s,t)}{\partial s^2}/\left(1+ \frac{\partial x(s,t)}{\partial s}\right)^{7/2}$$ The initial conditions at$t=0$are $$\text{init1}= x(s,t)\text{/.}t\to 0==0$$ $$\text{init2} = \frac{\partial x(s,t)}{\partial t}\text{/.}t\to 0==0$$ The boundary conditions at specified values of s are $$\text{bound1}= x(s,t)\text{/.}s\to 0==0$$ $$\text{bound2}=\frac{\partial ^2x(s,t)}{\partial t^2}\text{/.}s\to 1==-1 + 5\left/\left(1+ \frac{\partial x(s,t)}{\partial s}\text{/.}s\to 1\right)^{5/2}\right.$$ The condition bound2 can be simplified. In fact, defining the function $$\text{bound2a}=\ g(t)=\frac{\partial x(s,t)}{\partial t}\text{/.}s\to 1$$ which at the same time is the simplified boundary condition, bound2 reads an an ODE for$g$$$\text{eqg}=\frac{\ d g(t)}{\ d t}==-1 + 5\left/\left(1+ \ g(t) \right)^{5/2}\right.$$ The initial condition for$g\$ is taken from init2 to be

$$\text{initg}=g(0) == 0$$

Formulation in Mathematica

In order to have the boundary conditions complete we will solve first the ODE eqg. It is convenient to do this without imposing an initial condition. This will be done subsequently numerically. Here we go

solg = DSolve[g'[t] == -1 + 5/(1 + g[t])^(5/2), g[t], t]

(*
Out[150]= {{g[t] -> InverseFunction[
1/10 (-2 5^(2/5) Sqrt[10 - 2 Sqrt[5]]
ArcTan[(5 - 5 Sqrt[5] + 4 5^(4/5) Sqrt[1 + #1])/(
5 Sqrt[2 (5 + Sqrt[5])])] +
2 5^(2/5) Sqrt[2 (5 + Sqrt[5])]
ArcTan[(5 + 5 Sqrt[5] + 4 5^(4/5) Sqrt[1 + #1])/(
5 Sqrt[10 - 2 Sqrt[5]])] +
4 5^(2/5) Log[5 - 5^(4/5) Sqrt[1 + #1]] -
5^(2/5) (1 + Sqrt[5]) Log[
5 - 1/2 5^(4/5) (-1 + Sqrt[5]) Sqrt[1 + #1] + 5^(3/5) (1 + #1)] +
5^(2/5) (-1 + Sqrt[5]) Log[
5 + 1/2 5^(4/5) (1 + Sqrt[5]) Sqrt[1 + #1] + 5^(3/5) (1 + #1)] +
10 #1) &][-t + C[1]]}}
*)


In order to comply with initg we have to adjust the constant C[1] properly. This is done numerically thus

u0 = u /. FindRoot[0 == ((g[t] /. solg) /. t -> 0 /. C[1] -> u), {u, 1}]

(*
Out[213]= 1.19413
*)


The function is therefore (called gg henceforth)

gg[t_] := (g[t] /. solg /. C[1] -> u0)

Plot[gg[t], {t, -1, 5}, PlotRange -> {-1, 1},
PlotLabel -> "The function gg(t)"]
(* 150602_plot _gg.jpg *)


Now we can write down the complete equations and conditions

eq = D[x[s, t], {t, 2}] == D[x[s, t], {s, 2}] 1/(1 + D[x[s, t], s])^(7/2);

init1 = x[s, 0] == 0;

init2 = (D[x[s, t], t] /. t -> 0) == 0;

bound1 = x[0, t] == 0;

bound2a = (D[x[s, t], t] /. s -> 1) == gg [t];


And find the numerical solution

sol = NDSolve[{eq, init1, init2, bound1, bound2a}, {x[s, t]}, {s, 0, 1}, {t,
0, 1}]

(*
Out[224]= {{x[s, t] -> \!$$\* TagBox[ RowBox[{"InterpolatingFunction", "[", RowBox[{ RowBox[{"{", RowBox[{ RowBox[{"{", RowBox[{"0.", ",", "1."}], "}"}], ",", RowBox[{"{", RowBox[{"0.", ",", "1."}], "}"}]}], "}"}], ",", "\<\"<>\"\>"}], "]"}], False, Editable->False]$$[s, t]}}
*)


Here's the graph of x(s,t)

Plot3D[Evaluate[x[s, t] /. sol], {t, 0, 1}, {s, 0, 1}, PlotRange -> All,
PlotLabel -> "The solution x(s,t)", AxesLabel -> {"t", "s", "x(s,t)"}]
(* 150602_Plot3D _xst.jpg *)


Conclusion

The problem can be solved in Mathematica once we analyse carefully the equation and initial and boundary conditions (which is not trivial in this case). Despite the first impression based on a linearization the solution is not trivially zero. On the other hand is looks rather innocent and there is not much "suspense" in it.

• @ bbgodfrey : maybe you're right. But let's see, I'm just writing ... – Dr. Wolfgang Hintze Jun 2 '15 at 15:52
• In going from bound2 to eqg, you appear to have replaced $$\frac{\partial x(s,t)}{\partial s}$$ by $$\frac{\partial x(s,t)}{\partial t}$$. This can be verified by plotting the original boundary condition, Derivative[0, 2][x][1, t] - (-1 + 5/(1 + Derivative[1, 0][x][1, t])^(5/2)) using your numerical solution. – bbgodfrey Jun 2 '15 at 18:01
• @ bbgodfrey : you are right, I have made this wrong replacement. I'm trying to mend my solution. Let's see if I can impose bound2 directly ... No, as expected I get the message from NDSolve : Boundary condition bound2 should have derivatives of order lower than the differential order of the partial differential equation. So the problem of the OP seems to be ill defined. Continued in next comment. – Dr. Wolfgang Hintze Jun 2 '15 at 18:29
• @ bbgodfrey : A brute force modification, replacing the second derivative in bound2 by the first, leads, apart from a message of inconsitent boundary and initial conditions, to the trivial solution. So, the OP should be modified to a reasonable one by the author. Meanwhile I have shown a rather complicated way of writing down the number 0. ;-) – Dr. Wolfgang Hintze Jun 2 '15 at 18:31
• Yeah,I think so,This problem is not normal But fortunately the last condition could be reduced to a more simple one. :) – WateSoyan Jun 3 '15 at 2:00

Although the final boundary condition in the question certainly is unusual, the problem is solvable as follows by a do-it-yourself method of lines:

n = 100; h = 1/n; vars = Table[u[i, t], {i, 1, n}];
eqs = Join[
Table[D[u[i, t], {t, 2}] == (u[i - 1, t] - 2 u[i, t] + u[i + 1, t])/h^2
/(1 + (u[i + 1, t] - u[i - 1, t])/(2 h))^(7/2), {i, 1, n - 1}] /. u[0, t] -> 0,
{D[u[n, t], {t, 2}] == -1 + 5/(1 + (u[n, t] - u[n - 1, t])/h)^(5/2)}];
init1 = Table[u[i, 0] == 0, {i, 1, n}];
init2 = Table[(D[u[i, t], t] /. t -> 0) == 0, {i, 1, n}];
sol = NDSolveValue[{eqs, init1, init2}, vars, {t, 0, 1}];

ListPlot3D[Table[sol[[i]] /. t -> j h, {j, n}, {i, n}], DataRange -> {{0, 1}, {0, 1}},
PlotRange -> All, AxesLabel -> {"s", "t", "x"}, LabelStyle -> Directive[Bold, 12]]


Note: n = 100 provides better than 1% accuracy.

In answer to the comment below by Dr. Wolfgang Hintze, I have plotted our respective answers for s = 1 and 0 < t < 9:
The blue curve is his, and the brown one mine. They differ by about 25% at t = 1 and by much more at larger t. Both numerical solutions eventually fail, his at about t = 10, and mine at t = 36, although mine shows signs of numerical instability for t > 16. Presumably, better choices of integration parameters would extend both solutions.