13

This might seem extreme, but this can be analysed exactly. Take the system $$ \dot x_1 = x_2, \qquad \dot x_2=-x_1, \qquad x_1(0) = 1, \qquad x_2(0)=0. $$ Let $X=(x_1,x_2)$ be the state vector, $\delta t$ the time step, and $X^+$ the state vector for the next time step. Then the implicit Euler scheme is $$ X^+ = \delta t\left(\begin{array}{cc}0&1\\-1&...


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


11

If you want Matrix classes with an intuitive interface All the LAPACK and BLAS features Easy to learn and use API Easy to install Then I recommend you to have a look at my library FLENS. I designed it for exactly these kind of tasks. However, it requires a C++11 conform compiler (e.g. gcc 4.7 or clang). FLENS gives you exactly the same performance as ...


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

The standard approach is to use a self-starting time-marching algorithm with sufficiently small timestep (such that the order of accuracy is not spoiled) and compute the 5 non-initial value previous solutions. These are then used to "start" the BDF formula.


8

A simplification - the Crank-Nicolson method uses the average of the forward and backward Euler methods. The backward Euler method is implicit, so Crank-Nicolson, having this as one of its components, is also implicit. More accurately, this method is implicit because $u^{n+1}_i$ depends on $F^{n+1}_i$, not just $F^{n}_i$ This is means that the state at ...


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


6

The two properties are usually called positivity-preserving and monotonicity-preserving (makes it easier to find this question). Looking at http://www.ams.org/journals/mcom/2006-75-254/S0025-5718-05-01794-1/S0025-5718-05-01794-1.pdf and http://homepages.cwi.nl/~willem/DOCART/SIAM_HRS.pdf it seems that implicit Euler is the exception among implicit linear ...


6

I think you're mixing up ideas. The best thing here is to know the foundations, and it then "implicit FEM" will be a trivial idea (which is why there won't be a book specifically about that). Finite Element Methods are a form of spatial discretization where you use some kind of structured or unstructured mesh to re-write the problem as a system of equations ...


6

The reason is that GMRES can only be used for solving linear equations, i.e. equations of the form $Ax=b$, where $A$ is some matrix and $x,b$ are vectors. What GMRES does, essentially, is it approximates multiplication by the matrix $A^{-1}$ using a matrix polynomial of $A$. In this case (I assume) $f(y^{n+1},t)$ is not necessarily linear in the vector $y^{...


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

You would use backward euler method to solve a differential equation of the form $u_t =f(u,t)$ where $f$ is not necessarily a linear function in u. When f is non-linear, then the backward euler method results in a set of non-linear equations that need to be solved for each time step. Ergo, Newton-raphson can be used to solve it. For example, take $$\...


5

The implicit Euler method is unconditionally stable alright, but what you are doing is not the implicit Euler method. Rather, what you do is compute where the particle would be at the end of the time step using only 2-particle interactions, compute the location update, and then sum these updates up for all particle interactions. But the implicit Euler method ...


5

Your second one is the correct form. BDF approximates derivative with backward difference. Write $$ M(y) \dot{y} = f(y,t) $$ as $$ \dot{y} = M^{-1}(y) f(y,t) =: F(y,t) $$ write the BDF for this $$ \sum_{k=0}^s \alpha_k y_{n-k} = h F(y_n,t) = h M^{-1}(y_n) f(y_n,t_n) $$ Hence you get $$ M(y_n) \sum_{k=0}^s \alpha_k y_{n-k} = h f(y_n,t_n) $$


4

Backwards Euler, as you state you know, is time advances using the equation: $U_{t+\Delta t} = U_{t}+ U'(t,U_{t+\Delta t}) \Delta t$. You can not explicitly evaluate $U'_{t+\Delta t}$ since you don't know $U_{t+\Delta t}$. To make it more clear the equation is: $X = U(t) + U'(t,X)\Delta t$ Therefore you must solve the equation for X (which is U at the ...


4

This particular example is often called linearly implicit Euler. Its linear stability is identical to nonlinearly implicit Euler, but the nonlinear stability can be a limiting factor, especially for larger time steps. You can find some discussion for reaction-diffusion systems in Ropp, Shadid, and Ober 2004. More formally, this is the simplest example of a ...


4

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


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


4

To fix notation, denote the Butcher tableau by $$ \begin{array}{c|c} c & A\\ \hline & b^T \end{array} $$ where $b$ and $c$ are vectors of length $s$ (the number of stages) and $A$ is a $s \times s$ matrix. Consider the ODE $$ y' = f(t, y) $$ and suppose that $y_n$ is a given approximation to $y$ at time $t = t_n$. A general Runge-Kutta method (...


3

The Crank-Nicolson method is: $\frac{u^{n+1}_{i}-u^{n}_{i}}{dt} = \frac{1}{2}(F^{n+1}_{i}+F^{n}_{i})$ This method calculates the next state of the system, i.e. $u^{n+1}_{i}$, by solving an equation involving the previous states and the next state. In the case of the heat equation for example we would get a linear system and if we are using finite elements ...


3

Your question is correct, but this is still an open problem (see paper here). There are many methods to evolve CFL value in case of implicit methods. Every method has its own advantages and disadvantages. Choose the strategy which satisfies your needs or otherwise you can introduce efficient new strategy and publish it. All the best.


3

Define $$ \begin{align} A^{n+1} & = e\left[j^{n+1} - f\, θ_H^{n+1}\sinh\left(\frac{g\,n_A^{n+1}}{T}\right)\right] \\ B^{n+1} & = a\left[bP^{n+1}\,(1 - θ_H^{n+1} )^2 - c~(θ_H^{n+1})^2 - f\,θ_H^{n+1}\,\sinh\left(\frac{g\,n_A^{n+1}}{T}\right)\right] \\ C^{n+1} & = h(i - P^{n+1}) + d\,T\,(P^{n+1}-P_a)\left[bP^{n+1}~(1 - θ_H^{n+1} )^2 - c~(θ_H^{...


2

Basically the same question popped up on SO: What are the most widely used C++ vector/matrix math/linear algebra libraries, and their cost and benefit tradeoffs? (It adds some value to Geoff's answer.)


2

Since pdepe accepts systems of PDEs through vector-valued capacity, flux, and source terms, one way to accommodate your request would be to set the fluxes for all of the $\rho$ variables equal to zero. The capacity terms for the $\rho$ variables will all be 1, and the source terms for each variable are the non-flux terms on the right-hand side (the $\rho$ ...


2

As I said in my answer to your previous question, it's probably better if you don't try to write your own ODE solver, which is what you are doing now. There are a lot of good libraries out there that will solve a system of ordinary differential equations much better than if you roll your own implementation. If you can use them, you should; since your ...


2

It looks like the model you're trying to solve is: \begin{align} (1/\alpha(w,c))T_{t}(r,t) &= T_{rr}(r,t) + (p/r) \cdot T_{r}(r,t) \\ w_{t}(r,t) &= -(k_{1}(T(r,t)) + k_{2}(T(r,t)) + k_{3}(T(r,t)))w(r,t) \\ g_{t}(r,t) &= k_{1}(T(r,t))w(r,t) \\ a_{t}(r,t) &= k_{2}(T(r,t))w(r,t) \\ c_{t}(r,t) &= k_{3}(T(r,t))w(r,t) \end{align} where: $r =$...


2

One of the best ways to test a PDE solver is to use the method of manufactured solutions. Essentially, you modify the PDE (and discretization) by adding a source term that yields an exact solution known in advance. You can then compare your numerical solution against the exact solution for debugging purposes. You should test your numerical solution against a ...


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


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