12

Short answer If you only want second order accuracy and no embedded error estimation, chances are that you'll be happy with Strang splitting: half-step of reaction, full step of diffusion, half step of reaction. Long answer Reaction-diffusion, even with linear reaction, is famous for demonstrating splitting error. Indeed, it can be much worse, including "...


9

I assume, that you have conducted a space discretization, so that you are about solving the (vector-valued) ODE $$ \dot u_h(t) = F_h(t,u_h(t)), \text{ on [0,T] }, u_h(0) = \alpha. $$ via a numerical scheme $\Phi$ that advances the approximation $u_h^n$ at the current time instance $t=t^n$ to the next value $u_h^{n+1}$ at $t=t^{n+1}:=t^n+\tau$. Then your ...


8

Yes it is the time integration but it also means that: You have to solve a linear system of the type Ax=b in the implicit scheme where as in the explicit scheme you do not, as the lumped mass matrix only has diagonal entries so inv(M) is trivial. Your time step in the explicit scheme is limited by the CFL criteria for stability. Implicit schemes are ...


8

If you just slap together an implicit and an explicit method you will likely have order loss. You can do so with low order methods though, and Crank-Nicholson mixed with some other integrator is an easy way to get a decent second order integrator. Higher order IMEX integrators like the Kennedy and Carpenter Additive Runge-Kutta methods or the SBDF schemes ...


7

The FEM method for transient problems typically uses the method of lines, i.e. the spatial discretization is decoupled from the time discretization: \begin{equation} u^h(x,t) = \mathbf{\Phi}(x)^T \, \mathbf{U}(t) \end{equation} where $\mathbf{U}(t)$ is the vector of nodal quantities, assumed as unknown functions of time. Under this assumption the space-time ...


5

The terms "explicit" and "implicit" arise in the time discretization, and these terms are already used in the literature on ordinary differential equations (i.e., they are not specific to the finite element method). It would be worth taking a look at a book discussing the numerical solution of ODEs, e.g. Hairer & Wanner.


5

In general, just because method A provides certain guarantees (such as unconditional stability, energy conservation, being symplectic) does not imply that it is more accurate. In fact, a common observation is that the opposite may be true: For example, if a method is symplectic, then it guarantees that the error is zero with regard to certain quantities (e.g....


5

If we take your ODE, $$\frac{dy}{dx}=-\frac{x^2}{y},$$ multiply both sides by $y$ and integrate up, we see that the solutions look like $$ y^2 = C-\frac{2}{3} x^3. $$ Taking your initial condition, your real trajectory is then $$ y = \sqrt{1- \frac{2}{3} x}. $$ This is only real valued for the domain $-\infty < x<=1.5 $, so it's no real surprise that ...


4

You don't have just a first-order ODE so you cannot use an explicit Runge-Kutta method. Because of the square term, you cannot bring this into mass matrix form even. Instead, what you have is an implicit ODE. This falls into the class of problems known as Differential-Algebraic Equations, and that still works even though you don't have any pure algebraic ...


4

This is exactly the case when the lack of information in the question allows to answer it pretty certainly: it is certainly possible. The error would depend on many factors, including the conditioning of the original problem, particular details of the numerical implementation, and chosen simulation parameters. I do not see any contradiction yet. However, I ...


3

You have written down the midpoint rule as a two-step method, a member of the family of multi-step methods. For these methods, one can show that a multi-step method $$\alpha_{k} y_{n+k} + \alpha_{k-1} y_{n+k-1} + \ldots + \alpha_{0} y_{n} = h\left( \beta_{k} f_{n+k} + \ldots + \beta_{0} f_{n}\right)$$ is called stable if the polynomial $$\rho(z) = \alpha_{k} ...


3

In a one-step method $$ y_{n+1}=y_n+h\Phi_f(x_n,y_n,h) $$ one gets a truncation error for the exact solution $$ y(x_{n+1})=y(x_n)+h\Phi_f(x_n,y(x_n),h)+h^{p+1}\tau(x_n) $$ For the error propagation of $e_n=y_n-y(x_n)$ this gives approximately $$ e_{n+1}=e_n+h\partial_y\Phi_f(x_n,y(x_n),h)e_n-h^{p+1}\tau(x_n)+... $$ The error is then of order $p$, $e_n=c(x_n)...


2

Short Answer: There is no general result that would hold for all implicit schemes. The reason is that how your method behaves with respect to numerical diffusion depends on the specific combination of numerical time stepping scheme and spatial finite difference. However, note that if your discretisation is consistent with your PDE and your PDE does not have ...


2

You cannot merely adjust an explicit RK scheme into an implicit one, implicit routines are much more involved because each intermediate slope can depend upon slopes 'in the future'. This introduces a (possibly nonlinear) system of equations that needs to be solved within each RK step. That is nowhere to be found in the explicit method you posted above. If ...


2

All you need to do is find the stability criteria for whatever scheme (in this case Leap Frog Method) with respect to the simple model equation: $$\frac{dz}{dt} = \lambda z$$ for complex $\lambda$. Once you have the stability criteria, you need to find the eigenvalues of the matrix you have in your ODE system. Then just make sure each eigenvalue satisfies ...


2

I agree with Bill Barth that "time to solution" is a decent metric, but I'm afraid it has a couple of downsides: It's not portable: what's fast on one machine is not guaranteed to be fast on another (modern architectures differ not just in their floating point performance, but also in how their caches work or what their memory bandwidth is). It's ...


2

Backward Euler and the implicit trapezoidal rule are both unconditionally stable for this problem. If you're seeing instability then you haven't implemented them correctly.


1

Using the exact solution to estimate $\lambda$ may not give a time step that works for the numerical scheme. The local estimate of $\lambda$ is $$ \lambda = 400 t u $$ If you try this step $$ h = \min(2/20, 2/(400 t u+10^{-12})) $$ it is stable, but of course accuracy is poor. Here is a code import numpy as np import matplotlib.pyplot as plt t = 0.0 u = 1....


1

You seem to have given the 1D equations for the discretizations, even though the problem is in 2D. Regardless, the explicit method requires the least memory since you don't even have to form a matrix to compute the solution at the next time step. If you have the solution vector at time $t_i$, you can apply simple stencil operations and directly obtain the ...


1

In essence it does not matter what $\vec{j}$ means, but what it is. All worth to know is that your conservation equation is something like: $$\frac{\partial m}{\partial t} +\textrm{div}\,\vec{j}(m) = S(m) \tag{*}$$ The flux $\vec{j}(m)$ is a vector quantity that usually (and in this case it does) depends on the spatial difference in the variable $m$, $i.e.$ ...


1

Time to solution with (approximately) equal error is the best metric you have for comparing algorithms like this, but you need to do it across a range of mesh sizes and PDE coefficients. The explicit methods will have time step constraints that are $O(h^p)$ where $h$ is the mesh size and $p$ is related to the discretizations you use. This means that they ...


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