Problems with NDSolve and partial differential equations of several variables

Question

Suppose we have the following partial differential equation:
$$ 0 = \frac{ \partial w }{ \partial \tau } + \left( w + \sqrt{ h + \beta } \right) \frac{ \partial h }{ \partial \chi } $$ where $w$ = $w(\chi,\tau)$ and $h$ = $h(w)$ is implicitly given by: $$ w = 2 \left( \sqrt{ h + \beta } - \sqrt{ 1 + \beta } \right) + \sqrt{ \beta } \left( \ln \lvert \frac{ \sqrt{ h + \beta } - \sqrt{ \beta } }{ \sqrt{ h + \beta } + \sqrt{ \beta } } \rvert - \ln \lvert \frac{ \sqrt{ 1 + \beta } - \sqrt{ \beta } }{ \sqrt{ 1 + \beta } + \sqrt{ \beta } } \rvert \right) $$ Similar to the method of characteristics, we want to find paths in $(\chi,\tau)$ where $w$ = constant, which can be described by: $$ \frac{ d\chi }{ d\tau } = w + \sqrt{ h + \beta } $$ which allows us to define a general solution for the first equation given by: $$ w = \mathcal{F}\left( \chi - \left( w + \sqrt{ h + \beta } \right) \tau \right) $$ We also have the following boundary conditions and initial values (some helped to determine the first equation above):
$$ \begin{align} \lim_{\lvert x \rvert \rightarrow \infty} w\left( \chi, \tau \right) & = 0 \\ \lim_{\lvert x \rvert \rightarrow \infty} h\left( w\left( \chi, \tau \right) \right) & = 1 \\ \lim_{\tau \rightarrow 0} h\left( w\left( \chi, \tau \right) \right) & = 1 + \frac{ 1 }{ 2 } \cosh^{-1}{ \left( \frac{ \chi }{ 20 } \right) } \end{align} $$

I have tried several iterations on a theme, where I am trying to determine $w$ and $h$ as a function of $\chi$ and $\tau$, without success. The result I am looking for is the profile of the $w$ and $h$ versus position, $\chi$, and then be able to plot snapshots of that profile at different times, $\tau$.

Some variations resulted in the NDSolve::litarg error when I had originally tried to explicitly write $h$ = $h(w(\chi,\tau))$. So then I tried to reform the equations such that $w$ = $w(\chi,\tau)$ and $h$ = $h(\chi,\tau)$, but now I get NDSolve::underdet errors. I also tried using all of the equations (see code below), but then I get an over constrained system.

For instance, Mathematica was upset when I used 0 == Limit[w[x, t], x -> Infinity] instead of 0 == w[Infinity, t]. Unfortunately, I have run into a few errors, the source of which I did not fully understand, which resulted in my starting the typical I'm frustrated, so trial-and-error it is brute force method.

I assume I have an ignorance(inexperience)-based syntax issue, so any suggestions to help me see the error of my ways would be greatly appreciated.

Mathematica Code Trials
Mathematica version = 8.0.1.0
Operating System = Mac OS X x86 (v10.8.5)

hh[x_, t_] := h[w[x[t], t]];
ww[x_, t_] := w[x[t], t];

eq0 = ((w == -2 Sqrt[1 + bb] + 2 Sqrt[h + bb] - 
        Sqrt[bb] Log[(-Sqrt[bb] + Sqrt[1 + bb])/(
          Sqrt[bb] + Sqrt[1 + bb])] + 
        Sqrt[bb] Log[(-Sqrt[bb] + Sqrt[h + bb])/(
          Sqrt[bb] + Sqrt[h + bb])]) /. {h -> hh[x, t], 
      w -> ww[x, t]}) /. {h[w[x[t], t]] -> h[x, t], 
    w[x[t], t] -> w[x, t]};
eq1a = ((x[t] == (w + Sqrt[h + bb]) t) /. {h -> hh[x, t], 
      w -> ww[x, t]}) /. {h[w[x[t], t]] -> h[x, t], 
    w[x[t], t] -> w[x, t]};
eq1b = ((x'[t] == w + Sqrt[h + bb]) /. {h -> hh[x, t], 
      w -> ww[x, t]}) /. {h[w[x[t], t]] -> h[x, t], 
    w[x[t], t] -> w[x, t]};
eq2 = (ReleaseHold[
    Hold[(0 == D[w, t] + (w + Sqrt[h + bb]) D[w, x])] /. {h -> 
       h[x, t], w -> w[x, t]}]);
eq3 = (ReleaseHold[
    Hold[(0 == D[h, t] + w D[h, x] + h D[w, x])] /. {h -> h[x, t], 
      w -> w[x, t]}]);
eq4a = 0 == Limit[w[x, t], x -> Infinity];
eq4b = 0 == Limit[w[x, t], x -> -Infinity];
eq4c = 0 == Limit[w[x, t], Abs[x] -> -Infinity];
eq5a = 1 == Limit[h[x, t], x -> Infinity];
eq5b = 1 == Limit[h[x, t], x -> Infinity];
eq5c = (ReleaseHold[
      Hold[(h == 1 + 1/2 Sech[x/20])] /. {h -> hh[x, 0], 
        w -> ww[x, t]}] /. {x[0] -> x}) /. {h[w[x, 0]] -> h[x, 0]};
eq5d = ReleaseHold[Hold[(Derivative[1, 0][h][0, t] == 0)]];

xmax = 150.;
tmax = 25.;
soln = NDSolve[
  ({eq0, eq2, h[Infinity, t] == 1, h[-Infinity, t] == 1, 
     w[Infinity, t] == 0, w[-Infinity, t] == 0, eq5c, 
     eq5d} /. {bb -> 2.}),
  {w[x, t], h[x, t]},
  {t, 0, tmax},
  {x, -xmax, xmax},
  MaxStepSize -> 0.01, MaxSteps -> 10^6, SolveDelayed -> True
  ]  

[The folowing is only for house-keeping purposes, so ignore at your lesure]
Mathematica Code Trials: Things that didn't seem to work

(* 
eq5a=((ReleaseHold[Hold[(h==1)]/.{h->hh[x,t],w->ww[x,t]}]/.{x[t]->\
Infinity})/.{h[w[Infinity,t]]->h[0,t],h[w[-Infinity,t]]->h[0,t]})/.{h[\
w[x[t],t]]->h[x,t],w[x[t],t]->w[x,t]};
eq5b=((ReleaseHold[Hold[(h==1)]/.{h->hh[x,t],w->ww[x,t]}]/.{x[t]->-\
Infinity})/.{h[w[Infinity,t]]->h[0,t],h[w[-Infinity,t]]->h[0,t]})/.{h[\
w[x[t],t]]->h[x,t],w[x[t],t]->w[x,t]};
eq4a=((ReleaseHold[Hold[(w==0)]/.{h->hh[x,t],w->ww[x,t]}]/.{x[t]->\
Infinity})/.{h[w[Infinity,t]]->h[0,t],h[w[-Infinity,t]]->h[0,t]})/.{h[\
w[x[t],t]]->h[x,t],w[x[t],t]->w[x,t]};
eq4b=((ReleaseHold[Hold[(w==0)]/.{h->hh[x,t],w->ww[x,t]}]/.{x[t]->-\
Infinity})/.{h[w[Infinity,t]]->h[0,t],h[w[-Infinity,t]]->h[0,t]})/.{h[\
w[x[t],t]]->h[x,t],w[x[t],t]->w[x,t]};
eq2=(ReleaseHold[Hold[(0==D[w,t]+(w+Sqrt[h+bb])D[w,x])]/.{h->hh[x,t],\
w->ww[x,t]}])/.{h[w[x[t],t]]->h[x,t],w[x[t],t]->w[x,t]};
eq3=(ReleaseHold[Hold[(0==D[h,t]+w D[h,x]+h \
D[w,x])]/.{h->hh[x,t],w->ww[x,t]}])/.{h[w[x[t],t]]->h[x,t],w[x[t],t]->\
w[x,t]};
eq4a=ReleaseHold[Hold[(Limit[w,x->Infinity]==0)]/.{x->x[t],h->hh[x,t],\
w->ww[x,t]}];
eq4b=ReleaseHold[Hold[(Limit[w,x->-Infinity]==0)]/.{x->x[t],h->hh[x,t]\
,w->ww[x,t]}];
eq5a=ReleaseHold[Hold[(Limit[h,x->Infinity]==1)]/.{x->x[t],h->hh[x,t],\
w->ww[x,t]}];
eq5b=ReleaseHold[Hold[(Limit[h,x->-Infinity]==1)]/.{x->x[t],h->hh[x,t]\
,w->ww[x,t]}];
*)

(*
soln=NDSolve[
({eq0,h[Infinity,t]==1,h[-Infinity,t]==1,w[Infinity,t]==0,w[-Infinity,\
t]==0,eq5c,eq5d}/.{bb->2.}),
{w[x,t],h[x,t]},
{t,0,tmax},
{x,-xmax,xmax},
MaxStepSize->0.01,MaxSteps->10^6,SolveDelayed->True
];
soln=NDSolve[
({eq0,eq4a,eq5a,eq5c,eq5d}/.{bb->2.}),
{w[x,t],h[x,t]},
{t,0,tmax},
{x,-xmax,xmax},
MaxStepSize->0.01,MaxSteps->10^6,SolveDelayed->True
];
soln=NDSolve[
{eq0,eq1a,eq1b,eq2,eq3,eq4c,eq5a,eq5c,eq5d},
{w[x,t],h[x,t]},
{t,0,tmax},
{x,-xmax,xmax},
MaxStepSize->0.01,MaxSteps->10^6,SolveDelayed->True
];

soln = NDSolve[
{eq0, eq1a, eq1b, eq2, Limit[w[x[t], t], x -> Infinity] == 0, 
 Limit[w[x[t], t], x -> -Infinity] == 0, 
 Limit[h[w[x[t], t]], x -> Infinity] == 1, 
 Limit[h[w[x[t], t]], x -> -Infinity] == 1, 
 Limit[h[w[x[t], t]], t -> 0] == 1 + 1/2 Sech[x[t]/20], x[0] == 0},
{w[x[t], t], h[w[x[t], t]]},
{t, 0, tmax},
MaxStepSize -> 0.01, MaxSteps -> 10^6, SolveDelayed -> True
]
*)

Answer 1

This question relates to nonlinear wave steepening and I see now that I did not list everything that was necessary to address the problem. I am also quite certain that I phrased this question poorly on account of my novice-level experience with NDSolve in Mathematica.

The question spawned from my attempt to reproduce results from this paper. On the second page, they go into some details on how they analyzed the problem by simplifying their original quasi-general magnetized fluid equations. On the right-hand side of this page, the authors seemed to imply that they produced the results in their Figure 1 using the approach I outlined in my OP. Unfortunately, in my frustration I failed to catch an important detail on the last page of their paper.

The authors analyzed two equations numerically, shown here: $$ \begin{align} \left( \partial{\scriptstyle_{\tau}} + w \ \partial{\scriptstyle_{\chi}} \right) w + \partial{\scriptstyle_{\chi}} \left( h + \beta \ln h \right) & = 0 \\ \partial{\scriptstyle_{\tau}} \ h + h \ \partial{\scriptstyle_{\chi}} \ w + w \ \partial{\scriptstyle_{\chi}} \ h & = 0 \end{align} $$ where $h$ represents a normalized magnetic field component and $w$ a normalized velocity component. They assumed the following boundary conditions: $$ \begin{align} \lim_{\lvert x \rvert \rightarrow \infty} w\left( \chi, \tau \right) & = 0 \\ \lim_{\lvert x \rvert \rightarrow \infty} h\left( w\left( \chi, \tau \right) \right) & = 1 \\ \end{align} $$ and assumed the following initial conditions: $$ \lim_{\tau \rightarrow 0} h\left( w\left( \chi, \tau \right) \right) = 1 + \frac{ 1 }{ 2 } \cosh^{-1}{ \left( \frac{ \chi }{ 20 } \right) } $$ I finally realized, after wasting a considerable amount of time, that the pieces I was missing were: the first two PDEs in this answer and the initial condition for $w(\chi,0)$. The last missing piece is given by: $$ w = 2 \left( \sqrt{ h + \beta } - \sqrt{ 1 + \beta } \right) + \sqrt{ \beta } \left( \ln \lvert \frac{ \sqrt{ h + \beta } - \sqrt{ \beta } }{ \sqrt{ h + \beta } + \sqrt{ \beta } } \rvert - \ln \lvert \frac{ \sqrt{ 1 + \beta } - \sqrt{ \beta } }{ \sqrt{ 1 + \beta } + \sqrt{ \beta } } \rvert \right) $$ where we substitute our initial condition for $h(\chi,0)$ shown above into this equation. Then, we need only let Mathematica do the work. The code below works, though it is computationally expensive and made my computer work enough to get the fans going.

tmax = 100.;  (* Max. time step *)
xmax = 500.;  (* Max. spatial point *)
dtmaxf = 0.0001;  (* Max. fractional step *)
beta = 2.;  (* value for plasma beta *)
xfac = 20; (* normalization factor for Sech *)
xbounds = {-xmax/2, xmax};
tbounds = {-1, tmax};
(* Define initial condition functions *)
hatt0[x_, xo_] := 1 + 1/2 Sech[x/xo];
wofhb[hh_, bb_] := Sqrt[bb] + Sqrt[1 + bb] (-2 Sqrt[1 + bb] + 2 Sqrt[hh + bb] -
                   Sqrt[bb] Log[(-Sqrt[bb] + Sqrt[1 + bb])/()] + Sqrt[bb]
                   Log[(-Sqrt[bb] + Sqrt[hh + bb])/(Sqrt[bb] + Sqrt[hh + bb])]);
(* Define partial differential equations *)
eq0 = D[w[x, t], t] + w[x, t] D[w[x, t], x] + D[(h[x, t] + be Log[h[x, t]]), x] == 0;
eq1 = D[h[x, t], t] + w[x, t] D[h[x, t], x] + h[x, t] D[w[x, t], x] == 0;
eq2 = h[x, 0] == hatt0[x, xfac];
eq3 = w[x, 0] == (wofhb[h, be]) /. {h -> (hatt0[x, xfac])};
eq4a = h[xmax, t] == 1;
eq4b = w[xmax, t] == 0;
eqall = {eq0, eq1, eq2, eq3, eq4a, eq4b};

(* Find solutions to partial differential equations as function of [x,t] *)
soln = Block[{c = 0}, NDSolve[(eqall /. {be -> beta}),
             {w[x, t], h[x, t]},
             {x, xbounds[[1]], xbounds[[2]]},
             {t, tbounds[[1]], tbounds[[2]]},
             MaxSteps -> 10^6, SolveDelayed -> True, 
             Method -> {"MethodOfLines", "TemporalVariable" -> t}, 
             MaxStepFraction -> dtmaxf, SolveDelayed -> True, 
             StepMonitor :> If[Mod[c, 1000] == 0, Print["Step # = ", c]; c++, c++]]]

Then we can plot these results using the following:

(* Define time steps and colors for plots *)
tsolns = {0, 10, 20, 30, 40, 50};
cols = {Black, Purple, Blue, Green, Orange, Red};
(* Define plot outputs *)
plots = Table[Plot[((Evaluate[(h[x, t] /. Flatten[soln][[2]])]) /. {t -> tsolns[[i]]}), 
              {x, -100, 180}, PlotRange -> All, 
              PerformanceGoal -> "Quality", 
              PlotStyle -> {AbsoluteThickness[1], cols[[i]]}], 
              {i, 1, Length[tsolns], 1}];
(* Show plots *)
Show[plots, ImageSize -> {600, 600}, ImageMargins -> 0, AspectRatio -> Full]

The results are shown in the image below. Note that the last (red) time step shows a "fuzzy" part near x ~ 150. This is an error (perhaps not the right word) due to the steepening reaching a gradient catastrophe. At that point, the results become multivalued and I am not clever enough with Mathematica to know how to deal with this. I assume it has to do with the Stiffness settings in NDSolve, but again, I was not sure. I was just happy to get this working as well as it did. If any of the Mathematica whisperers on this site have a more efficient suggestion and a better way to handle the gradient catastrophe, I would greatly appreciate hearing about it.

Plot Description
The plot shows six snapshots of $h(\chi,\tau)$ at $\tau$ = 0 (black), 10 (purple), 20 (blue), 30 (green), 40 (orange), and 50 (red) versus $\chi$.

Update
After I updated java to Version 7 Update 71, the above code stopped working (i.e., it would quit the kernel and undefine everthing as if I had just opened the notebook) for Mathematica v8.0.1 on both my laptop (running OS X 10.9.5) and my desktop (running OS X 10.8.5). I updated my laptop to Mathematica 10 over the weekend and tried again. Then the Mathematica code worked again. My only guess as to the source of the problem was a change in the garbage collection of java, but JLink said that I was still using an earlier version of java. I wrote a formal report to Wolfram last Monday, but never heard back so I "bit the bullet" and upgraded hoping that would solve the problem. It did work, but I would have appreciated knowing this 5 days earlier so I would not have wasted a week I did not have trying to fix the issue.

I am not sure if the java update caused the problems, but I could not see why the bash update would have any effect on Mathematica. Then again, I am not an expert, so take my comments with a grain of salt.

Update 2
After upgrading to version 10, I discovered that Mathematica was trying to be too careful. Well, the support staff at Wolfram told me, but anyways the code needed to be altered to prevent a memory overflow. The simple solution was to do the following:

eqall = N[{eq0, eq1, eq2, eq3, eq4a, eq4b}];

Apparently, using what Mathematica defines as exact numbers caused a problem. Mathematica would finish the computation, but the memory would continue to disappear. After the addition of N[], the issue went away. Interestingly, at the point where the wave started to turn over, the code now fails. It had always ran into stiff solutions, as it should. But before it did not hit this limit until t ~ 51, whereas now it runs into problems at t > 47.95. In any case, I made a GIF that shows the steepening evolve, which is what I was going for in the first place.