Instability while NDSolving a wave equation

I'm trying to use NDSolve to solve a wave equations to check if it is easier and/or faster to use it instead of my old characteristics eq. method implementation.

I'm getting a lot of instability that I don't get with the characteristics method, and since these are simple equations, I wonder what is wrong... (hopefully, not the physical aspect of the problem...)

ans = Flatten@NDSolve[{
u[t, x]*D[d[t, x], x] + d[t, x]*D[u[t, x], x] + D[d[t, x], t] == 0,
D[d[t, x], x] + u[t, x]/9.8*D[u[t, x], x] + 
 1/9.8*D[u[t, x], t] + 0.0001 u[t, x]*Abs[u[t, x]] == 0,
u[0, x] == 0,
d[0, x] == 3 + x/1000*1,
u[t, 0] == 0,
u[t, 1000] == 0
},
d, {t, 0, 1000}, {x, 0, 1000}, DependentVariables -> {u, d}
]

Animate[Plot[(d /. ans)[t, x], {x, 0, 1000}, 
        PlotRange -> {{0, 1000}, {0, 6}}], {t, 0, 1000}
]

Can someone help me?

EDIT:

I've placed the NDSolve solution (following JxB's editing) with my characteristics solution, together on the same animation. They match close enough, with the exception of the initial quick oscillations. With time they tend do start do desynchronize, but I believe this is probably due to a small simplification we have to admit when deducing the characteristics.

Red: NDsolve; Blue: "manual" characteristics method;

press F5 (refresh your browser), to restart the animation from t=0.

(xx scale is the number of points I used with my "manual" method,开发者_高级运维 where each point represents 20 units of the NDSolve/physical scale)

Playing with NDSolve grid sampling, renders completely different oscillation effects. Does anyone have or know of a technique to ensure a proper integration?

By changing your coefficients to infinite precision (e.g., 1/9.8->10/98), and setting WorkingPrecision->5 (a value of 6 is too high), I no longer get the error message:

ans = Flatten@
  NDSolve[{D[u[t, x] d[t, x], x] + D[d[t, x], t] == 0, 
    D[d[t, x], x] + u[t, x] 10/98*D[u[t, x], x] + 
      10/98*D[u[t, x], t] + 1/10000 u[t, x]*Abs[u[t, x]] == 0, 
    u[0, x] == 0, d[0, x] == 3 + x/1000, u[t, 0] == 0, 
    u[t, 1000] == 0}, d, {t, 0, 1000}, {x, 0, 1000}, 
   DependentVariables -> {u, d}, WorkingPrecision -> 5]

Animate[
 Plot[(d /. ans)[t, x], {x, 0, 1000}, 
  PlotRange -> {{0, 1000}, {0, 6}}], {t, 0, 1000}]

I'm don't know this equation, so I don't believe the solution: small-scale oscillations grow initially, then are damped out.