30

In the solution of nonlinear hyperbolic PDEs, discontinuities ("shocks") appear even when the initial condition is smooth. In the presence of discontinuities, the notion of solution can only be defined in the weak sense. The numerical velocity of a shock depends on the correct Rankine-Hugoniot conditions being imposed, which in turn depends on numerically ...


17

The central question is which physical processes (waves or source terms) have time scales that you are interested in resolving and which you would prefer to step over. If you are not interested in the fastest time scale in the system, then the equations are called "stiff". Hyperbolic conservation laws are typically written as first-order systems $$ u_t + \...


14

The most common method is to reset negative values to some small, positive number. Of course, this is not a mathematically sound solution. A better general approach that may work and is easy, is to reduce the size of your time step. Negative values often arise in the solution of hyperbolic PDEs, because the appearance of shocks can lead to oscillations, ...


11

The principal numerical difficulty in solving a nonlinear first-order system of hyperbolic PDEs like the Euler equations (for compressible, inviscid flow) is that discontinuities (shock waves) appear in the solution after finite time, even if the initial data are smooth. In order to deal with this, most modern codes use both slope (or flux) limiters, which ...


11

If it's shock-capturing that you're interested in, I would suggest you use the finite volume method instead of the finite element method. When applied naively, FEM is actually notoriously bad at resolving shocks -- usually there are spurious oscillations or unwanted diffusion. Provided your original PDE is a conservation law, the FVM method will preserve ...


10

Here is a 97-line example of solving a simple multivariate PDE using finite difference methods, contributed by Prof. David Ketcheson, from the py4sci repository I maintain. For more complicated problems where you need to handle shocks or conservation in a finite-volume discretization, I recommend looking at pyclaw, a software package that I help develop. ""...


8

You could have a look at Fenics, which is a python/C framwork which allows quite general equations to be solved using a special markup language. It mostly uses finite elements though, but worth a look. The tutorial should give you an impression of how easy it can be to solve problems.


7

The short answer is: it requires specific work for different equations, but there are some general techniques that suggest how to do it. Essentially, given a first order evolution PDE $$u_t = Au + Bu$$ where $A,B$ are some (possibly differential) operators, the steady states are those for which $$Au + Bu=0.$$ It is common to use a splitting approach in ...


6

Assuming we are solving hyperbolic equations without any source terms and assuming we provide physical initial conditions, making sure the numerical scheme we use is Total Variation Diminishing is a good way of ensuring the "physicality" of the computed solution. Since a TVD scheme preserves monotonicity, no new minima or maxima will be created and the ...


6

I think one answer to your question is that certain communities simply always used conservative schemes and so it has become part of "the way it's done". One may argue whether that's the best way to do it, but that's about as fruitful as asking the British to drive on the right because it would simply be more convenient to have only on standard side. That ...


6

There is not much point using an implicit method for pure wave propagation because you have to resolve phase to have an accurate method. If you have a hyperbolic system in which some waves are very stiff (not interesting except for their influence on evolution of a slow manifold), you might want an implicit method. It is fairly problem-dependent whether you ...


6

There is definitely literature on schemes like this. Two keywords are Modified method of characteristics Semi-Lagrangian schemes After 20 minutes of googling: some possibly important papers are http://dx.doi.org/10.1137/0719063 and http://dx.doi.org/10.1137/0728024 (search forward from there). Those probably aren't the best references out there, but ...


6

A nonconvex flux function means that the characteristic velocity associated with a given characteristic field may not be a monotone function of the conserved variables. This can lead to non-classical Riemann solutions; for instance the total number of (shock, rarefaction, and contact) waves arising in the solution may be different from the number of ...


6

This is an important and challenging issue. Yes, using quadratic interpolation means that your solution values may lie outside the interval in which the initial data lie. This is not what we usually mean when we refer to numerical instability, but it is a potentially undesirable feature. Yes, forcing the interpolated values to lie in an interval destroys ...


6

Let me give a part of your answer, I would need some more indications from your side to answer you fully. So please read, and write some comments so that I can complete my answer. About notations in numerical methods There are a few mistakes in the way you write your equations. These are only details, but for someone who is used to it, it can a bit ...


6

Such problems (sometimes called lateral Cauchy problems) are in general not well-posed (meaning they either lack a solution, or there are infinitely many of them, or the solution is unstable under perturbations of the boundary conditions). For parabolic (or dissipative) equations, it makes sense to study the stationary limit (simply omit the term $u_t$ in ...


6

The stationary equation you show transports information from the right to the left via the advection term; it also diffuses slightly. If you switch off the diffusion term altogether, then you only have transport from the right to the left, and you need to also drop the boundary condition at the left: because information is from the right to the left, nothing ...


5

Godunov's method has an exact Riemann solver so no entropy fix is needed. A Roe solver (of which there are a few variants) uses a local linearization which has no diffusion to "fill in" the rarefaction fan, so it needs an entropy fix. Other approximate Riemann solvers, including Lax-Friedrichs, Rusanov, and the HLL family are inherently diffusive and do not ...


5

Starting from where David Ketcheson left me in his answer, a little bit more search revealed some historical notes. The scheme I outlined above was considered already back in 1900 by J. Massau, in Mémoire sur l'intégration graphique des équations aux dérivées partielles. The work is republished in 1952 by G. Delporte, Mons. The first (albeit brief) ...


5

As with many high order methods, the accuracy of the scheme is often less sensitive to the Riemann solver. None of the DG papers for hyperbolic problems will actually be using averages, however. The most common choice is a Rusanov (aka. Local Lax-Friedrichs) flux, which is very simple if you have an upper bound for the fastest wave speed.


5

I suggest looking at the literature on DG methods for incompressible flow, which has the mixed hyperbolic-elliptic character you mention. There are a lot of approaches. This paper, for instance, even uses an exact Riemann solver. This one suggests using a discontinuous space for the hyperbolic part and a continuous one for the elliptic part.


5

If the amplification factor is $\le 1$, then your scheme is stable; in your case it is exactly one, so, for this PDE, you won't have unbounded growth of errors. The Lax theorem then guarantees convergence of the method. In this narrow sense, the Crank-Nicholson method will work for you. (Have a look at Leveque's book if you want a good reference.) But, your ...


5

I have implemented the solution derived below using PyClaw in an IPython notebook. If you download that, you can adjust the initial values and see the computed solutions. General setup In the solution of a Riemann problem for the 1D Euler equations, there are generally 3 waves. Two of them are genuinely nonlinear; the other is a contact discontinuity, ...


5

Eigenvalues with zero eigenvalue correspond to purely oscillatory modes. You can see it by diagonalising the system. Your matrix $A$ can be written as \begin{equation} A = P \Sigma P^{-1} \end{equation} where $\Sigma$ is a diagonal matrix with entries $\lambda_n$ corresponding to the eigenvalues of $Q$. You can now transform your ODE to \begin{equation} dQ/...


5

First-order hyperbolic equations model conservation laws; as the alternative name "transport equations" suggests, they transport information along so-called "characteristic curves" with a finite speed of propagation. For the simple model equation $u_t + u_x = 0$, you can show (e.g., by separation of variables) that the solution should be of the form $$u(t,x)...


5

Your equation can be written in the following fashion (any spatial derivative approximation is valid), once space is discretised: $$\frac{1}{c}\frac{du_i}{dt}=-\left(\frac{\partial u}{\partial x}\right)_i(t) + v_i(t) \tag{*}$$ Keep in mind that $v_i(t) = v(x_i,t)$. Now the system of equations depends only on time $t$ you can apply Crack Nicholson method to ...


4

For the numerical solution of hyperbolic PDEs the use of Riemann solvers are essential components of conservative shock capturing methods for accurate simulation of wave problems which may have shocks (discontinuous jumps in conserved variables). To be able to obtain accurate solutions to such problems, we need to use proper upwinding techniques -- the ...


4

The above answers apply to time-dependent problems, but you could also demand positivity in a simple elliptic equation. In this case, you could formulate it as a variational inequality, giving bounds for the variables. In PETSc, there are two VI solvers. One uses a reduced-space method, where variables in active constraints are removed from the system to be ...


4

It is appears to be a standard finite difference discretization via method of lines of the wave equation $u_{tt} + c^2 \nabla^2u = 0, \quad (x,y)\in[0,L_x]\times[0,L_y], \quad t\geq0$ Standard finite difference approximations can be used to formulate a discrete problem by introducing the formally second-order accurate approximations $u_{tt}\approx \frac{u^...


4

PCHIP is not a conservative reconstruction, making it inappropriate for conservation laws. Furthermore, hyperbolic problems have discontinuous solutions so there is generally no benefit to a continuous reconstruction. Conservative monotone spline reconstructions are being investigated by the UK Met Office for use in tracer advection for atmosphere modeling, ...


Only top voted, non community-wiki answers of a minimum length are eligible