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


17

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


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


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

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

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

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

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

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


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

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

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


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

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

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

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

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

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

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

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


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