21

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


18

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


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

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

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


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


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


5

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


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


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

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

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

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


4

This is not an answer to your question, but more of an observation: More often than not, an ODE is "stiff" or a linear system is "nearly singular" because of a mistake either in deriving the equation to be solved, or in implementing it. Trying to find a way to solve what you have is then just a way to paper over the problem. If you were to find a way to ...


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


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

The number of rows and columns in the final global sparse stiffness matrix is equal to the number of nodes in your mesh (for linear elements). 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 ...


2

I assume you're trying to solve an equation that looks like: \begin{align} -\nabla \cdot (a(x)\nabla{u}) = f, \end{align} for $x$ in some domain $\Omega$, although the same approach would be fine (for a residual evaluation, anyway) if $a$ were also a function of $u$. The stiffness matrix will take the form \begin{align} A_{ij} = \int_{\Omega}a(x)\nabla\...


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

Odespy needs from you a 1st order system of ODEs, namely something of the form \begin{align} \frac{d\mathbf{y}}{dt} = \mathbf{f}(\mathbf{y},t), \tag{1} \end{align} where, for your particular problem, $y_i$ are single-variate complex-valued functions that depend solely on time (that means you need to arrive at the semi-discrete form first by discretizing the ...


2

I don't know of a better way to handle this with SciPy since I don't think it has event handling. But if you're willing to venture beyond SciPy, the following software have the capability, either documented as event handling or rootfinding: Sundials' CVODE MATLAB's ode23 and ode15s Julia's DifferentialEquations.jl Sundials wrapper and Rosenbrock methods ...


2

Rosenbrock methods utilize embedded lower order methods in order to calculate errors for adaptive time stepping. In addition, Rosenbrock methods do not have to solve an implicit system (just a linear system). There is no form of iteration then that takes place in them (unless you're using a Krylov linear solver). Maximum number of steps for stiff solver can ...


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


2

Due to your formulation, I call $X(z) = \begin{bmatrix} x(z) \\ p_{x}(z) \end{bmatrix}$ so your ODE is written in matrix form: $$X^{'}(z) = C(z) X(z)$$ Where: $C(z) = \begin{bmatrix} 0 & A(z) \\ B(z) & 0 \end{bmatrix}$. Your general formula by using backward Euler method is: $$\frac{X(z+\Delta z) - X(z)}{\Delta z} = C(z+\Delta z) X(z+\Delta z)$$ ...


2

See comments for the original discussion. The lack of L-stability of the trapezoidal rule in the first step is the source of your problem. A simpler and famous toy problem that shows the point is the Curtiss-Hirschfelder equation $$y'=-2000(y-\cos(t)) \\y(0)=0$$ Integrating this with Backward Euler and the Trapezoidal rule gives the following result (I wrote ...


1

The term "stiff" has many definitions, but a commonly understood one would be "there are two or more timescales associated with the process being modeled, and these timescales are very different". In your case, you have a single ODE whose solution is a single function $y(t)$. As a consequence, it has only a single timescale, and so one would not consider ...


1

If you want to stick to the Python scientific family, you could opt for Assimulo. Assimulo is a wrapper around a lot of ODE integrators, providing a common interface. If you happen to be running on Windows, you can find wheels at Christoph Gohlke's page. Assimulo allows to handle state events to allow for discontinuities in the RHS function of your ODE but ...


1

Mathematically speaking, stiffness is meaningless for a single differential equation, and is rather attributed to a set of differential equations that have different time-scales (e.g. when trying to solve two coupled equations with time-scales of 1 second and 1 day, respectively). However, a single equation can also be referred to as stiff if certain ...


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