23

Causality indicates that information only flows forward in time, and algorithms should be designed to exploit this fact. Time stepping schemes do this, whereas global-in-time spectral methods or other ideas do not. The question is of course why everyone insists on exploiting this fact -- but that's easy to understand: if your spatial problem already has a ...


21

J.M's comment is right: you can find an interpolating polynomial and differentiate it. There are other ways of deriving such formulas; typically, they all lead to solving a van der Monde system for the coefficients. This approach is problematic when the finite difference stencil includes a large number of points, because the Vandermonde matrices become ill-...


15

This is a well-framed question and a very useful thing to understand. Korrok is correct to refer you to von Neumann analysis and LeVeque's book. I can add a bit more to that. I'd like to write a detailed answer, but at the moment I only have time for a short one: With $\alpha=\beta=1/2$, you get a method that is absolutely stable for arbitrarily large ...


10

Generally speaking, you'll want to use an implicit method for parabolic equations (the diffusion part) -- explicit schemes for parabolic PDE need to have a very short timestep to be stable. Conversely, for the hyperbolic part (advection) you'll want an explicit method as it's cheaper and doesn't disrupt the symmetry of the linear system you have to solve by ...


10

The reason is that with the exception of linear problems, if you do a Fourier (or other) decomposition in time, you end up with a significant number of problems that are coupled globally in time. In other words, you have to solve lots of problems on the entire time interval concurrently. That will typically bust your computational or memory budget. The ...


9

Just as an aside, your github documentation is fantastic. This is just a guess from DG methods, which can have similar issues if numerical fluxes aren't chosen carefully (I figure FV methods are a subset of DG methods). If you're using interpolation from cell centers to define your fluxes, then this should be equivalent to using the average as a numerical ...


8

Similar to the causality Wolfgang mentioned in his post, we could see the reason why time dimension is special from Minkowski spacetime point of view:$ \newcommand{\rd}{\mathrm{d}} $ The $(3+1)$-dimensional spacetime has an inner product defined as $$ (A,B) = A_x B_x + A_y B_y + A_z B_z - \dfrac{1}{c^2}A_t B_t $$ if $A$ and $B$ are two 1-form in Minkowski ...


8

Many of us in scientific computing simply have well-equipped laptops for regular software development tasks, some multicore workstations for smaller-scale testing, and access to clusters for larger runs. To give you an idea: My laptop is a Dell M3800 (4-core Intel i7, hyperthreading, 16GB of RAM). This is good enough to regularly compile my software and do ...


7

The Galerkin approach (seeking an approximation from a given subspace $U$ such that the residual is orthogonal to another given subspace $V$) is indeed very general (and not restricted to finite-dimensional spaces). In the context of the numerical solution of partial differential equations, there are essentially two conditions that $U$ and $V$ have to ...


7

For simplicity, assume that there is only one parameter $t$ rather than your two. In order that you can have continuous eigenspaces, you need to assume that the associated eigenvalues do not nearly cross. (For nearly crossing eigenvalues it may very well happen that the eigenspaces are essentially exchanged though the eigenvalue curves do not touch. This ...


7

I would not insist on demanding a geometrical meaning from Galerkin methods in general. There is a connection, but it becomes less meaningful as you extend it further and further. (In a sense, it would be better to call Galerkin methods "generalized projection methods".) To really understand the connection between collocation and Galerkin methods requires ...


6

In the case you show, the solution has a boundary layer. If you can't resolve it because your mesh is too coarse, then for all practical matters the solution is discontinuous to the numerical scheme. Now, if you just apply a standard discretization to this problem, the discrete solution is the result of applying a linear projection operator to the exact ...


6

If the right hand side were independent of $u$ then one would generally use the averaged form $$ (1-\theta)s_(x,t^{n+1}) + \theta s(x,t^n). $$ In the nonlinear case you can't do that easily, as you note, but you can at least use some kind of extrapolation, for example approximate $$ (1-\theta)s(x,t^{n+1},u^{n+1}) + \theta s(x,t^n,u^n) \approx (1-\...


6

While there are some exceptions (e.g. fully discrete finite element methods), temporal discretization generally implies an inherently sequential dependence in flow of information. This dependency restricts semi-discrete algorithms (BVP in space, IVP in time) to compute solutions to subproblems in sequential manner. This discretization is usually preferred ...


6

I am currently implementing a VoF method (a geometrical method for two phase flow simulation on Eulerian mesh) that is native to structured grid, on an unstructured grid, so here are my experiences so far (please note that what I'm writing comes from working with a specific implementation): unstructured mesh: pros fast generation of meshes for complex ...


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

If you need Jacobian matrix information for a numerical method, you should calculate the Jacobian matrix of the discretized form of the equations, since that will be consistent with the discretized equations you are solving.


6

Every major class of discretization is "open-ended" in the sense that there are decisions with no obviously/provably correct answer in the general case, so some decisions are made based on how they perform for the target problem. Additionally, each major class has active research on new extensions. AMR has more choices than static-grid discretizations, so ...


6

I believe the best approach here is to use a threshold based on the local average brightness of the image. Setting the threshold to be 90% of the mean value of the 11x11 grid surrounding each pixel gives results that are about as good as you can expect with such a low resolution image. For each pixel you just need to compute the mean brightness of the ...


6

We have the following problem: $$\frac{\partial u}{\partial t}+v\color{red}{\frac{\partial u}{\partial x}}-\nu\color{blue}{\frac{\partial^2u}{\partial x^2}}=0 \tag{*}$$ The function $u$ may represent for example the concentration that propagates at velocity $v>0$ and disperses in a medium with viscosity $\nu>0$. Since only we are discussing how terms ...


6

Instead of directly integrating over the area, it is often more convenient to use the divergence theorem to replace the area integral with an integral over the boundary edges. The divergence theorem in three dimensions for a vector function $F$ can be written $$ \int_V \nabla \cdot {\bf F} dv = \int_S {\bf F} \cdot ds $$ That is, an integral over the ...


5

Others have said it all, but I just wanted to add a simple, yet sometimes subtle, point. Your upwind discretization remains conservative as long as you use a consistent interpolation of $a(x)$ on the cell boundaries. What I mean by consistent is that the only condition that the interpolation needs to satisfy is $$ a_{i+1/2}^+ = a_{i+1/2}^- $$ In other ...


5

Adaptive refinement. There are even optimal error estimators for your exact problem, though you can't go too wrong on this problem by dividing cells where the magnitude of the gradient is largest.


5

Section 2.5 of this PDF document goes through some additional details and the error does in fact work out to be second order. The key is that the term multiplying the $\Delta y^2/\Delta x$ term is equal to $$\Delta x\left( \dfrac{\partial^4u}{\partial x\partial y^3} +\mathcal{O}(\Delta x^2) \right)$$ The $1/\Delta x$ cancels the $\Delta x$ in this term ...


5

Full space-time discretization of time-dependent partial differential equations is indeed a thing. If you use a structured mesh in time (in the sense that the time discretization does not depend on space) and appropriate choice of trial and test functions, you can fit several standard time-stepping methods (Crank-Nicolson, implicit Euler or some Runge-Kutta ...


5

If you wish to know only the change in density you are interested in the dilatation or volumetric strain of the stress tensor associated with the transformation. If you also wish to know about the anisotropy of the stretching you should look at computing the right Cauchy-Green deformation tensor, $\Delta$, and it's associated eigenvalues and eigenvectors. ...


5

For the particular equation you are solving (called the minimal surface equation), the functional you are trying to minimize is $$ J(u) = \int_\Omega \sqrt{1 + |\nabla u|^2} \; dx. $$ You can find a derivation of the equations, as well as a discussion of solution approaches in lectures 31.5 and following here: http://www.math.tamu.edu/~bangerth/videos....


5

The concept of $p$-nonconforming meshes can also apply to continuous FEM, not just DG. For continuous FEM, continuity is enforced by enforcing a single set of degrees of freedom on a face shared between two elements. If those elements have varying polynomial degrees, the trace space on the face must be made the same. This may be done by restricting the ...


5

You can't. Nonlinear systems of equations are in general not solvable exactly. What you need to do is to use a method to solve nonlinear systems, of which there are of course quite a lot: A simple approach would be to use $D(U)\approx D(U^{n-1})$, where $D^{n-1}$ is the solution of the previous time step. A possibly smarter approach would be to use $D(U)\...


4

MATLAB and Octave are susceptible to the same subtle floating-point issues that Python is, where you can get a slightly unexpected result if you do not anticipate rounding issues. On the other hand, this is quite convenient syntax to have! You can easily "bake" a simple function to do this from the existing functionality in numpy.linspace, which is a very ...


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