146

I've gathered the following from online research so far: I've used Armadillo a little bit, and found the interface to be intuitive enough, and it was easy to locate binary packages for Ubuntu (and I'm assuming other Linux distros). I haven't compiled it from source, but my hope is that it wouldn't be too difficult. It meets most of my design criteria, and ...


42

Newton's method may not converge for many reasons, here are some of the most common. The Jacobian is wrong (or correct in sequential but not in parallel). The linear system is not solved or is not solved accurately enough. The Jacobian system has a singularity that the linear solver is not handling. There is a bug in the function evaluation routine. The ...


24

This document was written in March 2009 to help in the choice of a linear algebra library for a scientific library. It evaluates portability, high-level interface and licensing for several libraries, among them Eigen, GSL, Lapack++ MTL, PETSc, Trilinos and uBlas. It seems to be particularly fond of Flens and Seldon. (One of the requirements was that C++ ...


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


17

Of all the projects listed above, there are really only two heavy-weights that are extremely widely used (and for good reasons): PETSc and Trilinos. Both are professionally developed and have a large developer base. All the others are rather small projects compared to these two, and I would recommend going with them because (i) they will be supported for a ...


15

It's quite common in computational fluid dynamics to use implicit schemes similar to what you propose. The ones I know of are based on compact finite difference formulas (not simply on replacing $n$ with $n+1$ in existing schemes). For instance, one of the most widely used schemes was developed by Lele in 1992 in this paper with >2500 citations. Such ...


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


12

There is no reason that you cannot do what you wrote. One of the reasons that this is uncommon is that there for hyperbolic (advection) type problems the domain of dependence is finite. Thus an explicit methods makes sense from a computational efficiency standpoint. The implicit scheme you have written will require solving a linear system, albeit in the ...


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

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


6

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


6

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

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

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

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

Since this is a very active area of research, I hesitate to attempt an answer to this but I have some experience on what not to try. Do Not: Take old application code A and old application B, then try to couple them together Use archaic (implying unusable in hindsight) code, instead of building a new application Require a huge framework (> 10 required ...


5

I strongly advocate a fully coupled assembly since this can easily reproduce the operator split versions. Specifically, the routines which calculate the residual and Jacobian for different physics can be separate, but the framework should be able to combine them to form a unified residual for the entire system. This is how PETSc works. Then, the operator ...


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

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


5

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

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

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) $$


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