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 ...


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

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

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

You are correct: If you satisfy the CFL condition, then all that guarantees is that your scheme is stable, i.e., the numerical solution does not go to infinity. But the CFL condition says nothing about how accurate the numerical solution is. For that, indeed $\Delta z$ and $\Delta t$ must also be small enough compared to the features of the exact solution. ...


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

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

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

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

The CFL number (or simply Courant number) is defined locally: $\nu = u(x,t) \frac{\Delta t}{\Delta x}$. The velocity may vary in time and space and the grid may be non-uniform (in $x$ and/or $t$). So the CFL number can be different at each point in time and space. A necessary condition for convergence of a consistent method is that the CFL condition be ...


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

This is not an answer. It is an attempt to point out some of the errors in the question. Any source of FDM schemes for the transport equation starts by explaining that explicit central differences for the equation of the type $u_t+au_x=0$ can cause oscillations and unconditionally unstable. The explicit centered difference scheme for this equation (in ...


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, ...


4

In 1-D, the first discretization you've presented is correct. The matrix equation does not look right, though. For starters, it's not clear what your boundary conditions would be or how you would incorporate them. In $N$ dimensions, the advection equation looks like \begin{align} \frac{\partial{u}}{\partial{t}} + \sum_{i=1}^{N}\frac{\partial}{\partial{x^{i}...


4

The first two can be decomposed into a system of first-order advection equations of which there are a number of standard test problems often involving both smooth and discontinuous initial conditions and all of which can be solved analytically. For instance this solution is to the linearized acoustics equations. The choice between one or the other is ...


4

Central differencing schemes are not stable if you have advection dominated problems. There really was no other trivial [1] alternative to developing upwind schemes. [1] There are a few other stabilization methods, of course, but in the finite difference of finite volume context, almost everything that was developed over the first 30 years of numerical ...


4

To elaborate on Wolfgang's answer: since hyperbolic PDE semi-discretizations with centered differences have purely imaginary eigenvalues, they are only neutrally stable. For linear problems (e.g., the acoustic wave equation or linear Maxwell's equations) this is fine and such methods are commonly used. For instance, in electromagnetics centered differences ...


4

This is a Hamilton-Jacobi equation. You can read about how to apply WENO to such equations in Section 4 of Chi-Wang Shu's 2009 WENO review paper, and references therein.


4

I have two extra points I would like to add to Wolfgang's answer. A formulation of the CFL condition that I find more useful than the classic formula is this: A necessary condition for the stability of a numerical scheme is that the numerical domain of dependence bounds the physical domain of dependence. This is exactly what good old $$ \dfrac{\Delta ...


4

Your numerical solution is probably just getting more accurate as you increase the number of grid points. Do you know or have you tried to derive the analytic (exact) solution for this problem? By looking at your plots, it seems like the exact solution has a shock (discontinuity) occurring at around x = 4.5, and the numerical method is resolving it with ...


4

This doesn't answer your question directly but instead suggests an alternate strategy. In general, unless you are interested in experimenting with numerical methods, I recommend using an existing PDE solver instead of trying to write one yourself. This advice is particularly true for numerically-challenging problems, like yours appears to be. Since you are ...


3

I have not studied this particular system before, and I'm sure that someone who has could say much more than I will. I don't think there is any reason to expect that discretizing the equations in the form you have them will lead to highly accurate conservation of total energy. In this form, conservation of total energy depends on cancellation of certain ...


3

$\def\rmC{{\mathrm C}}$ $\def\vr{{\vec r}}$ $\def\vA{{\vec A}}$ $\def\ve{{\vec e}}$ $\def\l{\left}\def\r{\right}$ In the following I try to show that Bale's form is correct. If you complain, leave a comment. The equations we can trust are the integral equations. Mass Balance: $$ \l[\int_0^x \rho(x',t') A(x') dx'\r]_{t'=0}^t + \int_{t'=0}^t \l[(\rho u)(x'...


3

There are plenty of higher order schemes out there. But as per Godunov's theorem, only first order scheme can be monotone and hence not create oscillations. This resource gives a brief idea about the construction and analysis of finite difference schemes. In REA (Reconstruct, Evolve , Average) algorithm, required order polynomial is reconstructed and ...


3

Many times, the equations to be solved represent a physical conservation law. For example, the Euler equations for fluid dynamics are representations of conservation of mass, momentum, and energy. Given that the underlying reality that we are modeling is conservative, it is advantageous to choose methods that are also conservative You can also see ...


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