19

So there is a ton to say about this, and we will actually be putting a paper out that tries to summarize it a bit, but let me narrow it down to something that can be put into a quick StackOverflow post. I will make one statement really early and keep repeating it: you cannot untangle the efficiency of a method from the efficiency of a software. The details ...


12

The property follows from the property of the corresponding (weak form of the) partial differential equation; this is one of the advantages of finite element methods compared to, e.g., finite difference methods. To see that, first recall that the finite element method starts from the weak form of the Poisson equation (I'm assuming Dirichlet boundary ...


10

Stiffness involves some separation of scales. In general, if you are interested in phase of the fastest mode in the system, then you have to resolve it and the system is not stiff. But frequently, you are interested in the long-term dynamics of a "slow manifold" rather than the precise rate at which a solution off the slow manifold approaches it. Chemical ...


10

Part 1 Small eigenvalues are not included in the definition of stiffness for ODE (initial value problem) systems. There is no satisfying definition of stiffness that I know of, but the best definitions I've come across are: If a numerical method with a finite region of absolute stability, applied to a system with any initial conditions, is forced to use ...


9

You seem to be limited by stability, which is expected since diffusion is stiff as you refine the grid. Good methods for stiff systems are at least partly implicit. It will take some effort, but you can implement a simple multigrid algorithm (or use a library) to solve this system with a cost of less than ten "work units" (essentially the cost of one of your ...


9

A better way to look at it is that for a stiff problem, any stable explicit calculation leads to an error that is much smaller than the required error tolerance. There are many good methods for automatically detecting stiffness using explicit schemes, especially embedded Runge-Kutta pairs. See for example: Detecting stiffness with the Fehlberg (4, 5) ...


9

1. Can we numerically detect stiffness just by applying explicit methods? Suppose you have an initial value problem for some ODE on $[0,10]$. You take considerably large stepsize $\tau=1$ and an explicit Euler method, make your calculations with constant step size $\tau$ and get these points: You estimate the error and it appears to be big. Okay, then you ...


7

You will want to read up on operator splitting methods. In essence, in every "macro time step" you would treat fast processes by doing many "micro time steps" in one half of the algorithm, and then do a single macro time step for the slow processes in the other half. For higher order, you will want to use what's known as "Strang splitting".


6

To answer your questions: As far as I know, in practice, if explicit methods require extraordinarily small time steps relative to your time scale of interest (see answers to this question on what it means for an ODE to be stiff) in order to yield accurate results, then for all intents and purposes, your problem is stiff. To determine requirements on step ...


6

There are so many Runge Kutta methods, including Dormand-Prince 45 Cash-Karp 54 Fehlberge (sic) 78 Is there any comparison between them? Well, sure. Here are some traits to compare: Is the method implicit or explicit? (All of your examples are explicit RK methods.) What is the order of convergence? Are there any embedded error ...


6

The term that you want to search for is multiple timestepping (see, for instance, [1-3]). [1] http://www.cs.unc.edu/Research/nbody/pubs/external/Berne/tuckerman-berne-rossi91.pdf [2] http://www3.nd.edu/~izaguirr/papers/newM3paper.pdf [3] http://arxiv.org/abs/1307.1167


6

The only documentation I know about for the implementation of ode23t is in the paper which documents the implementation of ode23tb, the TRBDF2 method in MATLAB. As usually implemented, the trapezoidal rule is not strictly a one-step method because the truncation error estimate makes use of two previously computed solution values. The method is not ...


5

ode15s is designed to handle stiff systems of ODEs so I doubt if the problem you are encountering is that your "equations are too stiff" It is more likely that your spatial discretization has an error for some reason or your have some other MATLAB programming error. I suggest the following approach to debug this: Set the final integration time for ode15s to ...


5

There is no solid definition for stiff equations. I like Shampine's working definition the best: a differential equation is stiff if explicit methods are less computationally efficient than implicit methods. All of the other definitions about eigenvalues are purposefully vague because they cannot capture the fact that the "eigenvalues which matter" is ...


4

You can see the formulas for explicit RK4 at Mathworld. Given the simplicity of the expressions it is very reasonable to code one for yourself which is why you may not have found a module. However, if as you say you are solving a stiff IVP then an explicit Runge-Kutta method is not appropriate (the regions of absolute stability are bounded). Perhaps you ...


4

The zero eigenvalue is not particularly harmful. Think of the linear system $$ \dot x = A x, x(0) = x_0 $$ then the components of $x(t)$ decay in proportion to the (magnitude of the) eigenvalues of $A$ (assuming the eigenvalues are all non-negative), whereas the components of $x(t)$ in direction of eigenvectors that belong to zero eigenvalues simply stay ...


4

There are a couple of questions implicit in your post: How does one deal with non-uniqueness of the algebraic equations generated by any implicit numerical method? Typically you have a very good initial guess -- the previous time step solution. This will usually ensure that e.g. Newton's method converges to the correct root. For extremely large values ...


4

The approach of using two equations for the real and imaginary parts of an equation is often used. It may lead to more cumbersome formulas, but it is definitely possible and common. An example where this is done in deal.II (a library that I maintain) is here: https://www.dealii.org/developer/doxygen/deal.II/step_29.html


4

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


4

In BDF schemes for $\dot y = f$, one uses $$ f(t_n)=\dot y(t_n) $$ and tries to approximate $\dot y(t_n)\approx \sum_{j=0}^k\alpha_k y_{n-j}$ by the current value $y_n$ (that is to be computed) and the $k$ previously computed approximations. In the presented approach, in $(5)$, $y$ is approximated as a polynomial $p$ in $t$ fitted to $y_{n-j}$, so that the ...


3

The number of rows and columns in the final global sparse stiffness matrix is equal to the number of nodes in your mesh. For example if your mesh looked like: then each local stiffness matrix would be 3-by-3. Once all 4 local stiffness matrices are assembled into the global matrix we would have a 6-by-6 global matrix. For example the local stiffness matrix ...


3

There is the case of production and destruction of radioactive isotopes in nuclear fuel. You typically use nuclear fuel for more than 3 years (hence, you fission isotopes into radioactive daughter products that decay with certain half-lives that can go from microseconds to years or more). So during the irradiation, if you want to track all these isotopes (to ...


3

A differential equation $y'=F(y)$ may be called stiff in $[0,T]$ if, for some $y$, the matrix $TF'(y)$ has some huge eigenvalues with negative real part. Thus the simplest example is $y'=-y$. Indeed, the equation $y'=-10^6 y$ for $t\in[0,1]$ is obviously stiff, but $y'=-y$ is as stiff for $t\in[0,10^6]$. An explicit method such as Euler's needs the same ...


3

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


3

First of all, I cannot judge the quality of the DotNumerics implementation of Radau5, I can only assume it is a direct port of the Fortran original from Hairer and Wanner. A full description of their Radau5 implementation is in their book on solving stiff differential equations. The test you refer to is a check to see if the Jacobian should be recomputed. ...


3

Can we expect radau5 to cope with discontinuities, or should we integrate the two trajectories separately? Differential equation methods cannot easily handle discontinuities like this. If you step over a discontinuity, you cannot prove that you will not have order loss. In fact, you normally will have loss of accuracy. Because of this, you want to make sure ...


2

From a practical point of view: the A5 processor is not that much powerful, so you can wait a few HW iterations, or if your ipod / ipad are going to be connected to the internet, solve your problem remotely or in the cloud.


2

In fact Jed Brown has cleared the question for me. What I'm doing now is just putting his words in the context. Both 2d linear ODE systems from above are stiff (i. e. hard to solve with explicit methods) on relative big time intervals (e. g. [0,1]). The linear systems with large stiffness ratio can be considered "more stiff" because most likely one needs ...


2

If stiffness of element is not positive, then the system is not stable. So the model is most likely not correct. Look at the most basic equation of harmonic oscillator $$m x''(t) + k x(t) = f(t)$$ The solution is unstable if $k$ is negative (look at the roots of the characteristic equation). It means the solution will blow up. The stiffness has to be a ...


2

If you relax the independence between $\epsilon$ and $\Delta t$ to be "approximate", i.e. use a method which has a really large stable region even if it's not A-stable, then there are plenty of methods which can get good performance. A standard set of methods for this are the Backward Differentiation Formulas (BDF) or the Numerical Differentiation Formulas (...


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