21

van der Houwen's statement is correct, but it is not a statement about all fifth-order Runge-Kutta methods. The "Taylor polynomials" he is referring to are (as you seem to know) just the polynomials of degree $p$ that approximate $\exp(z)$ to order $p$: $$P_p(z) = \sum_{j=1}^p \frac{z^j}{j!}$$ For the fifth-order polynomial, it turns out that $|P_5(i\...


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


15

There are two main classes of solutions to be discussed in this regard. "Sufficiently" Smooth Solutions In Strang's classical paper it is shown that the Lax equivalence theorem (i.e., the idea that consistency plus stability implies convergence) extends to nonlinear PDE solutions if they have a certain number of continuous derivatives. Note that that ...


13

The CFL condition states that the "mathematical domain of dependence" must be (asymptotically) contained in the numerical domain of dependence. For hyperbolic problems, this provides a bound $\Delta t < C \Delta x$ that is useful at all resolutions. For a parabolic problem, it merely requires that $\Delta t \in o(\Delta x)$ in the limit $\Delta x \to 0$. ...


13

If the solution of $Ax=b$ is unstable, the matrix is very ill-conditioned (i.e., has a very large condition number), and (paraphrasing Lanczos) no amount of mathematical trickery can make it stable. The best you can hope for is to solve a different problem that is a) stable and b) gives you a solution that is sufficiently close; this is called regularization....


13

The singular value decomposition for a symmetric matrix $A=A^{T}$ is one and the same as its canonical eigendecomposition (i.e. with an orthonormal matrix-of-eigenvectors), while the same thing for a nonsymmetric matrix $M=U \Sigma V^T$ is just the canonical eigenvalue decomposition for the symmetric matrix $$ H=\begin{bmatrix}0 & M\\ M^{T} & 0 \end{...


12

The formula $$\mathrm{logsum}(x,y)=\max(x,y)+\mathrm{log1p}(\exp(-\operatorname{abs}(x-y))$$ should be numerically stable. It generalizes to a numerically stable computation of $$\log \sum_i e^{x_i} = \xi+ \log\sum_i e^{x_i-\xi},~~~\xi=\max_i x_i$$ In case the logsum is very close to zero and you want high relative accuracy, you can probably use $$\...


12

MAC finite difference schemes (Harlow and Welch 1965) are uniformly stable, but require smooth structured grids and are only second order accurate. Finite element methods are preferred for unstructured and high order methods. For continuous Galerkin finite element methods, there are no known spaces that have optimal approximation properties and are ...


11

The best approach is to use an ODE solver that is guaranteed to conserve the norm of the initial condition, i.e., for which $\|y_n\| = \|y_0\|$ for all $n\in\mathbb{N}$. Such solvers exist, and are called geometric integrators, since they preserve geometric properties of the exact solution (in this case, that energy is conserved, i.e., $\frac{d}{dt}\|y(t)\| =...


11

Stability does indeed mean that small changes in the data lead to small changes in the solution. This can be shown for the linear advection equation through the energy method; in proving the stability of a discretization, we seek to mimic the energy estimate of the exact problem. Considering the PDE, $$ \frac{\partial u}{\partial t} + a \frac{\partial u}{\...


10

Very short answer: for a comprehensive reference, you can't beat Hairer and Wanner volume II. Short answer: Here are some MATLAB scripts to plot the stability region of a linear multistep or Runge-Kutta method, given the coefficients. You could also use the Python package nodepy (disclaimer: it's my package and it's not the most polished piece of software, ...


9

You can just take Clenshaw's recurrence $$ u_k(x) = 2xu_{k+1}(x)-u_{k+2}(x)+\color{red}{a_k},\\ f(x) = x u_1(x)-u_2(x)+\color{red}{a_0} $$ and differentiate it directly: $$ u_k'(x) = 2xu_{k+1}'(x)-u_{k+2}'(x) + \color{blue}{2u_{k+1}(x)},\\ f'(x) = x u_1'(x)-u_2'(x) + \color{blue}{u_1(x)}. $$ Note that now the derivatives of partial sums, $u_k'(x)$ satisfy ...


8

It all boils down to finding the largest or smallest eigenvalue of the stiffness matrix. You can do that analytically for a uniform mesh, and you then arrive at the traditional formula $\tau \le ch^2$. For non-uniform meshes, the relationship between mesh size and eigenvalue is not obvious any more; in particular, there is no analytic formula for the ...


8

These are all standard questions discussed in most books on ODE solvers. I would recommend Hairer & Wanner.


8

You can apply linear stability analysis. That is, for given $u=(x,v)$ compute the linearization $Df(u)$ of the right hand side if the equation is $u‘=f(u)$. The problem is stiff if those differ by orders of magnitude. At a glance, I would not expect this. You can determine a good step size by running the problem again with half the size. If the results are ...


8

For a hyperbolic system of equations, you can write your equation as $$ \frac{\partial \mathbf{u}}{\partial t} + [\mathbf{A}] \frac{\partial \mathbf{u}}{\partial x} = 0$$ and then perform an eigendecomposition $\mathbf{A} = \mathbf{Q} \mathbf{\Lambda} \mathbf{Q}^{-1} $ where $ \mathbf{\Lambda} $ is a diagonal matrix of the eigenvalues, then defining $\...


8

Let's denote by $\otimes,\oplus,\ominus$ (I was lazy trying to get circled version of division operator) the floating-point analogs of exact multiplication ($\times$), addition ($+$), and subtraction ($-$), respectively. We'll assume (IEEE-754) that for all of them $$ [x\oplus y]=(x+ y)(1+\delta_\oplus),\quad |\delta_\oplus|\le\epsilon_\mathrm{mach}, $$ ...


8

Notice that $\hat C_*=(1-r F(h))\hat C_n$, but the sign in front of $r$ is lost when you use $\hat C_*$ inside $\hat C_{n+1}$. I'm looking at p.149 of Numerical Solution of Time-Dependent Advection-Diffusion-Reaction equations by Hundsdorfer and Verwer on Google books. (I'm going to use their signs.) The stability region of Heun's method is $$|g(z)| = \...


8

Numerical solution of the advection equation with centered differences in space and forward Euler in time is unconditionally unstable. So the behavior you are seeing is expected. Here is a nice lecture by Gil Strang ( MIT 18.086) where he discusses the instability of this method. He also shows how the simple centered difference method can be modified to ...


7

If you substitute, at least for your analysis, $\frac{\partial u}{\partial x}$ by $u_x$, you can write your system as $$ \begin{bmatrix} 0 & 0 \\ I & I \end{bmatrix} \frac{d}{dt} \begin{bmatrix} p_h(t) \\ u_{x,h}(t) \end{bmatrix}+\begin{bmatrix} -\partial_h & \partial_h \\ -\Delta_h & 0 \end{bmatrix}\begin{bmatrix} p_h(t) \\ u_{x,h}(t) \end{...


7

What you are looking for is what is called "Automatic error analysis" and is the subject of Chapter 26 of Higham's book "Accuracy and Stability of Numerical Algorithms", 2nd ed., SIAM Publishers. One technique he describes is using direct search optimization: try to formulate your problem as an optimisation problem and use the optimization algorithm to ...


6

I don't know about the theoretical complexity. But the practically fastest algorithm for finding all zeros of a polynomial given as a finite power series is the algorithm by Traub [1]. It factors a degree $n$ polynomial in typically $O(n^2)$ time. See also [2]. For factoring a polynomial accessible only through its function values, see, e.g., my paper [3]....


6

For the compressible Navier-Stokes equations, a convenient possibility is to choose the time step by considering the stability criterion of the inviscid (hyperbolic, Euler) and the viscous (parabolic, diffusive) terms independently from one another: For the Euler term you have "something like" $ \Delta t_\mathrm{Hyperbolic} \leq \displaystyle\min_{\mathbf{...


6

This is an important and challenging issue. Yes, using quadratic interpolation means that your solution values may lie outside the interval in which the initial data lie. This is not what we usually mean when we refer to numerical instability, but it is a potentially undesirable feature. Yes, forcing the interpolated values to lie in an interval destroys ...


6

I'm going to answer a more general question than the one you asked: do the eigenvalues of an initial value ODE determine the stability of the solution? Here I'm referring to mathematical stability, not numerical stability. Of course, a "yes" to this question is a necessary condition for a "yes" to your question. And unfortunately, the answer is "no". In ...


6

You're learning that the two forms (with $-\Delta u$ and with $-2\nabla \cdot \varepsilon(u)$) are not equivalent. They lead to different boundary terms in the bilinear form after integration by parts. In the first case, you get a term involving $n\cdot \nabla u$, in the latter $2n \cdot \varepsilon(u)$, and the "natural" boundary conditions you can enforce ...


6

How accurately? Is there a reason you're not just using numpy.mean? Is that not sufficient? If you need to compensate for floating-point error, you can try using Kahan summation. Here's some pseudocode (PDF). For a review of floating-point summation techniques see: N.J. Higham, SIAM, 1993.


6

From your link, consider the definition of the round off error and the statement "Since the exact solution must satisfy the discretized equation exactly, the error must also satisfy the discretized equation." This is actually only true if the PDE is homogeneous, that is, if we can write it in the form $\mathcal{L}(u;u_t,u_x,u_{xx},\ldots)=0$, with all terms ...


6

There are a few ways you can do this. Even though RK4 fails unless stepsizes are really small, one thing that can happen a lot is that the model can start out more stiff than it really is in the full timespan. Thus I would still try something adaptive like MATLAB's ode45 or Julia's Tsit5() before ruling out non-stiff solvers. If that fails, then I would try ...


6

But clearly, this is not the case as my programs do come up with (an approximate) solution though. I believe you did not continue the integration until you see that your integration is not convergent and is not bounded. I could rewrite your system of ODEs as: $$\dot{x_{1}} = x_{2}$$ $$\dot{x_{2}} = -kx_{1}$$ Or in matrix form: $$\dot{X} = AX$$ Where: ...


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