NDSolve: how to solve "ndnum error"

Question

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}]     

Answer 1

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 = NDSolve`FiniteDifferenceDerivative[Derivative[2], xgrid, DifferenceOrder -> 2]["DifferentiationMatrix"]; 
  M1 = NDSolve`FiniteDifferenceDerivative[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 = NDSolve`FiniteDifferenceDerivative[Derivative[2], xgrid, 
   DifferenceOrder -> 4]@"DifferentiationMatrix"; M1 = 
 NDSolve`FiniteDifferenceDerivative[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}}]]