# My discretization of a wave equation in first-order form does not give correct solutions. What should I do?

I haven't much experience with conservation laws, shocks, etc. After reformulating my wave equation to 1st order system (velocity-stress): $$\frac{\partial v}{\partial t} + A \frac{\partial v}{\partial x} + B \frac{\partial v}{\partial y} = 0$$ where: $v=(v_{x},v_{y},\sigma_{x},\sigma_{y},\sigma_{xy}) \\ A = A(\rho, \mu, \lambda)\\ B = B(\rho, \mu, \lambda)$

I put my Dirichlet BC to the first component (some constant vx) and obtained a stiff behavior. I've tried the same models with 2nd order equation and saw beautiful waves propagating. Here I get immediate propagation and loads of instabilities with explicit solvers (FD, RK4). The finer the mesh, the more high-frequency oscillations I get. Maybe I'm used to displacements and don't understand this formulation...

My question is: Is there any way of solving this kind of problem with ordinary tools or should I bother with Riemann solvers, etc. ? My long-term plan is to extend it to DG. There's plenty of papers on Godunov scheme, etc. but I couldn't quite catch what's going on though...

EDIT: Yes, it's a linear hyperbolic system. I used the term "stiff" in a more 'colloquial' manner. I can't find a stable explicit timestepping scheme for this system with Dirichlet BC. Even for extremely low dt (1e-20), way below CFL I get immediate waves in the whole domain. On the other hand, for implicit, the solution does not propagate correctly. For instance, sinusoidal excitation u1 = y = sin(100*t) with 2nd order system looks like this (mesh is obviously too coarse but quite fine): For 1st order system cosinusoidal excitation v1 = y' = 100*cos(100*t) like this: The same mesh size, the same timestepping (Crank-Nicolson). I suppose this is not physically correct...

EDIT: The wave does not propagate from the middle too: • According to what you have written, $A$ and $B$ do not depend on $\nu$. In that case, this is a 2D linear hyperbolic system (in non-conservative form). That problem isn't ordinarily stiff, so you'll need to edit the question and explain why you think it is stiff (or what you mean by "obtained a stiff behavior). – David Ketcheson Feb 12 '16 at 10:44
• One is missing information to guess what is the problem, e.g. if you discretize the spatial derivatives with central differences, you might observe such behavior. I would try to solve a problem with only one initial function being nonzero in a small finite support in the middle of domain to exclude the boundary conditions as the current source of your problems. For some time the boundary values should not influence the numerical solution. – Peter Frolkovič Feb 12 '16 at 14:39
• @bubupankowski - then I would guess the space discretization is one of the candidates for a source of your problems. It seems you fix now the value at boundary below another constant value in the initial state and your "standard discretization" that does not use information on the direction of wave propagation, creates the oscillations. Try to fix the boundary at the same constant value as the initial constant value and change the initial function only in the middle of domain (far away of boundary) to some small wave to see for short time how it develops. It might give you some insight. – Peter Frolkovič Feb 12 '16 at 15:28
• If you tried centered differences + RK4 and you don't see stable behavior even for small CFL numbers, then there is an error in the implementation. It's not really possible to say more based on the information given. – David Ketcheson Feb 12 '16 at 16:14
• @bubupankowski- in your last remark I am missing the mass matrix before v_old, i.e. in the final step instead of M * v_new = v_old + ... , it should be M * v_new = M * v_old+ ... – Peter Frolkovič Feb 12 '16 at 17:10