# NDSolve: how to solve "ndnum error"

I have a problem of finding the shape of a current. please see the details here NDSolve does not accept the boundary condition

I wrote a code but it gives the "ndnum" error. I tested the recommendation in the document for solving this error but it was not helpful. Also, I tried different values for t0 (initial condition) but still no results.

Any help or hint would be appreciated.

ClearAll[h, x, n, t]
hx = D[h[x, t], x];
ht = D[h[x, t], t];
n = 1;     (* put  1 for the start, later we will  test different n*)

sqrtTerm = Sqrt[1 + 4*n*Abs[hx]];
pde = Sign[-hx]*1/2*D[h[x, t]*(sqrtTerm - 1), x] == -ht

(*based on 2017 paper Zeng, we decided to put h[x, t0]=1*)
solution =
NDSolve[{Sign[-hx]*1/2*D[h[x, t]*(sqrtTerm - 1), x] == -ht,
h[x, 0.001] == 1, h[1, t] == 0, Derivative[1, 0][h][0, t] == 0},
h, {x, 0, 1}, {t, 0.001, 10}]

• Dear @Alex Trounev , I saw one of your answers to a question here mathematica.stackexchange.com/a/229198/96497 I am aware that my problem should be solved by the method of line. would you please take a look at my problem and give me some suggestion?
– AWer
Commented Feb 12 at 10:55
• Could you give a link to the 2017 paper Zeng? Commented Feb 12 at 12:00
• journals.aps.org/prfluids/abstract/10.1103/… If you do not have access I can send it to you by email
– AWer
Commented Feb 12 at 12:13
• Thank you. Is it original equation that you proposed or it taken from some paper? Commented Feb 12 at 12:40
• the formula is original, but it is in the context of Gravity currents. usually with the help of similarity solutions and scaling the pde is transferred to an ode. but here I need to solve the pde numerically and then compare it with the results of the scaled ode
– AWer
Commented Feb 12 at 12:46

We can solve this problem using the Euler wavelets collocation method and method of lines. Supposed that h[x,t] is a real function, we express Abs[hx] as Sqrt[hx^2], so we have

ClearAll[h, x, n, t]
hx = D[h[x, t], x];
ht = D[h[x, t], t];
n = 1;

sqrtTerm = Sqrt[1 + 4*n*Sqrt[hx^2]];
pde = Sign[-hx]*1/2*D[h[x, t]*(sqrtTerm - 1), x] == -ht;


Nevertheless it can't be solved with NDSolve due to some reason. Therefore we use wavelets as a method to discretize the equation on x

OEm[m_, x_] :=
Sqrt[2  m +
1] Sum[(-1)^(m - k) x^k  Binomial[m, k] Binomial[m + k, k], {k, 0,
m}]; UE[m_, t_] := OEm[m, t];
psi[k_, n_, m_, t_] :=
Piecewise[{{2^((k - 1)/2)  UE[m, 2^(k - 1)  t - n + 1], (n - 1)/
2^(k - 1) <= t < n/2^(k - 1)}, {0, True}}];
PsiE[k_, M_, t_] :=
Flatten[Table[psi[k, n, m, t], {n, 1, 2^(k - 1)}, {m, 0, M - 1}]];
k0 = 2; M0 = 4; With[{k = k0, M = M0},
nn = Length[Flatten[Table[1, {n, 1, 2^(k - 1)}, {m, 0, M - 1}]]]];
dx = 1/(nn); xl = Table[l*dx, {l, 0, nn}]; xcol =
Table[(xl[[l - 1]] + xl[[l]])/2, {l, 2, nn + 1}]; Psijk =
With[{k = k0, M = M0}, PsiE[k, M, t1]]; Int1 =
With[{k = k0, M = M0}, Integrate[PsiE[k, M, t1], t1]];
Int2 = Integrate[Int1, t1];
Psi[y_] := Psijk /. t1 -> y; int1[y_] := Int1 /. t1 -> y;
int2[y_] := Int2 /. t1 -> y;
wA = Table[wa[i][t], {i, nn}]; wB = Table[wb[i][t], {i, 2}];
w2[x_] := wA . Psi[x]; w1[x_] := wA . int1[x] + wB[[1]];
w0[x_] := wA . int2[x] + wB[[1]]  x + wB[[2]];
eqw = With[{w =
w0[x]}, -D[w, t] == (-(1/2))*
Sign[Derivative[1, 0][h][x,
t]]*(Derivative[1, 0][h][x,
t]*(-1 +
Sqrt[1 + 4*Sqrt[Derivative[1, 0][h][x, t]^2]]) + (2*
h[x, t]*Sign[Derivative[1, 0][h][x, t]]*
Derivative[2, 0][h][x, t])/
Sqrt[1 + 4*Sqrt[Derivative[1, 0][h][x, t]^2]]) /. {h[x, t] ->
w0[x], Derivative[1, 0][h][x, t] -> w1[x],
Derivative[2, 0][h][x, t] -> w2[x]}];
eqnw = Table[eqw, {x, xcol}];
icx = With[{w = w0[x]}, w == 1 /. t -> 10^-3];
ic = Table[icx, {x, xcol}];
bc = With[{w = w0[x]},
Join[{w1[x] == 0} /. x -> 0, {w == 0} /. x -> 1]]; varAll =
Join[wA, wB];
icn = Join[ic, bc /. t -> 10^-3]; eqn = Join[eqnw, D[bc, t]]; var1 =
D[varAll, t];
{vec, mat} = CoefficientArrays[eqn, var1];
f = Inverse[mat // N] . (-vec); vr0 = varAll /. t -> 10^-3;
{v0, mat0} = CoefficientArrays[icn, vr0];
sol0 = LinearSolve[mat0, -v0];
icn0 = Table[vr0[[i]] == sol0[[i]], {i, Length[vr0]}];
sol1 = NDSolve[{Table[var1[[i]] == f[[i]], {i, Length[var1]}], icn0},
varAll, {t, 10^-3, 10}];


Visualization

ParametricPlot3D[{t, x, w0[x] /. sol1[[1]]}, {t, 10^-3, 10}, {x, 0,
1}, ColorFunction -> Hue, AxesLabel -> {"t", "x", "h"},
Boxed -> False, BoxRatios -> {1, 1, 1}, PlotPoints -> {40, 35},
MaxRecursion -> 1, MeshStyle -> {Red, Blue}]


Update 1. This problem can be solved with FDM as follows

L = 1;  dx = L/20; xgrid = Range[0, L, dx]; nn = Length[xgrid];
M2 = NDSolveFiniteDifferenceDerivative[Derivative[2], xgrid, DifferenceOrder -> 2]["DifferentiationMatrix"];
M1 = NDSolveFiniteDifferenceDerivative[Derivative[1], xgrid, DifferenceOrder -> 2]["DifferentiationMatrix"];
wA = Table[wa[i][t], {i, nn}];
w1 = M1 . wA; w2 = M2 . wA;
rhs = (-(1/2))*Sign[Derivative[1, 0][h][x, t]]*Derivative[1, 0][h][x, t]*(-1 + Sqrt[1 + 4*Sqrt[Derivative[1, 0][h][x, t]^2]]) -
(h[x, t]*Derivative[2, 0][h][x, t])/Sqrt[1 + 4*Sqrt[Derivative[1, 0][h][x, t]^2]] /. {h[x, t] -> wA, Derivative[1, 0][h][x, t] -> M1 . wA,
Derivative[2, 0][h][x, t] -> M2 . wA};
eq = Table[D[wA[[i]], t] == -rhs[[i]], {i, 2, nn - 1}];
bc = {(M1 . wA)[[1]] == 0, wA[[-1]] == 0};
ic = Table[wA[[i]] == 1 /. t -> 10^(-3), {i, 2, nn - 1}]; icn = Join[ic, bc /. t -> 10^(-3)];
eqn = Join[eq, D[bc, t]]; var1 = D[wA, t];
{vec, mat} = CoefficientArrays[eqn, var1];
ff = Inverse[N[mat]] . (-vec); vr0 = wA /. t -> 10^(-3);
{v0, mat0} = CoefficientArrays[icn, vr0];
sol0 = LinearSolve[mat0, -v0];
icn0 = Table[vr0[[i]] == sol0[[i]], {i, Length[vr0]}];

Off[General::partd];
f[t_, x_] := Evaluate[ff /. Table[wa[i][t] -> x[[i]], {i, nn}]]

rk2[f_, h_][{t_, x_}] := Module[{k1, k2}, k1 = f[t, x];
k2 = f[t + h/2, x + h k1/2];
{t + h, x + h k2}];
tf = 2; dt = 1/1000; sol =
NestList[rk2[f, dt], {0, icn0[[All, 2]]}, Round[tf/dt]];


Visualization

h = Interpolation[
Flatten[Table[{{sol[[j, 1]], xgrid[[i]]}, sol[[j, 2, i]]}, {j,
Length[sol]}, {i, nn}], 1], InterpolationOrder -> 1]

Plot3D[h[t, x], {t, 0, 2}, {x, 0, 1}, ColorFunction -> Hue,
PlotTheme -> "Scientific", AxesLabel -> {"t", "x", "h"},
PlotRange -> All, PlotPoints -> 50, MeshStyle -> {Red, Blue}]


Note, in this example we used the Runge-Kutta second order algorithm on $$0\le t \le 2$$ with step $$10^{-3}$$ and difference scheme of second order on $$0\le x\le 1$$ with step 1/20.

Update 2 This is more advanced code with the Runge-Kutta 4th order algorithm on $$0\le t\le 10$$ and difference scheme of 4th order on $$0\le x\le1$$

 L = 1; tend = 10; dx = L/20; xgrid = Range[0, L, dx]; nn = Length[xgrid];
M2 = NDSolveFiniteDifferenceDerivative[Derivative[2], xgrid,
DifferenceOrder -> 4]@"DifferentiationMatrix"; M1 =
NDSolveFiniteDifferenceDerivative[Derivative[1], xgrid,
DifferenceOrder -> 4]@"DifferentiationMatrix";
wA = Table[wa[i][t], {i, nn}];
w1 = M1 . wA; w2 = M2 . wA;

rhs = (-(1/2))*Sign[Derivative[1,0][h][x,t]]*(Derivative[1,0][h][x,t]*(-1+Sqrt[1+4*Sqrt[Derivative[1,0][h][x,t]^2]])+(2*h[x,t]*Sign[Derivative[1,0][h][x,t]]*Derivative[2,0][h][x,t])/Sqrt[1+4*Sqrt[Derivative[1,0][h][x,t]^2]]) /. {h[x, t] -> wA, Derivative[1, 0][h][x, t] -> M1 . wA,
Derivative[2, 0][h][x, t] -> M2 . wA};

eq = Table[D[wA[[i]], t] == -rhs[[i]], {i, 2, nn - 1}];
bc = {(M1 . wA)[[1]] == 0, wA[[-1]] == 0};
ic = Table[wA[[i]] == 1 /. t -> 10^-3, {i, 2, nn - 1}]; icn =
Join[ic, bc /. t -> 10^-3];
eqn = Join[eq, D[bc, t]]; var1 = D[wA, t];
{vec, mat} = CoefficientArrays[eqn, var1];
ff = Inverse[mat // N] . (-vec); vr0 = wA /. t -> 10^-3;
{v0, mat0} = CoefficientArrays[icn, vr0];
sol0 = LinearSolve[mat0, -v0];
icn0 = Table[vr0[[i]] == sol0[[i]], {i, Length[vr0]}];

Off[General::partd];
f[t_, x_] := Evaluate[ff /. Table[wa[i][t] -> x[[i]], {i, nn}]]

rk4[f_, h_][{t_, y_}] := Module[{k1, k2, k3, k4}, k1 = f[t, y];
k2 = f[t + h/2, y + h k1/2];
k3 = f[t + h/2, y + h k2/2];
k4 = f[t + h, y + h k3];
{t + h, y + h/6*(k1 + 2 k2 + 2 k3 + k4)}]


Example of usage

dt1 = 1/1000; tf = 10; sol1 =
NestList[rk4[f, dt1], {0, icn0[[All, 2]]}, Round[tf/dt1]];


Visualization

h = Interpolation[
Flatten[Table[{{sol1[[j, 1]], xgrid[[i]]}, sol1[[j, 2, i]]}, {j,
Length[sol1]}, {i, nn}], 1], InterpolationOrder -> 8];
Plot[Table[h[t, x], {t, {.1, .5, 1, 2, 5, 10}}] // Evaluate, {x, 0,
1}, Frame -> True, FrameLabel -> {"x", "H"},
PlotLegends -> Table[Row[{"t =", t}], {t, {.1, .5, 1, 2, 5, 10}}]]


• This is amazing. thank you very much. I appreciate it. it will takes time for me to understand it as I am new in Mathematica. But just one question, As I had no idea about Euler wavelet method ,I searched a little bit and only found that its applicable in wave and signal problems. my formula is for the fluid propagates into a porous medium. what is the logic to use this method for such a problem? it would be grate if you explain a bit more.
– AWer
Commented Feb 12 at 15:38
• @AWer This kind of wavelets nothing to do with signal process. This is the system of orthogonal functions to implement the decomposition of solution. In this example we use only 10 collocation points to compute numerical solution. We can compare it with other methods like FEM an FDM. Commented Feb 12 at 17:23
• thank you very much
– AWer
Commented Feb 13 at 8:54
• Thank you very much. it is great. Would you please guide me how to get the 2d plot of the results? I need to have the profile of the current. H vs x plot that shows the evolution of the current in different times
– AWer
Commented Feb 15 at 10:22
• @AWer Could you specified the different times like for example times={0,1,2,3,4}? Commented Feb 15 at 11:51