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


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

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


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


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}{\...


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

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


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


7

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


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

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


5

On your first question: I assume that by "usual discretization matrix" you mean either the 3-point finite difference discretization in 1d, or what you get using linear finite elements. In either case, it's not actually the Crank-Nicolson scheme that determines this. It's true that for the two spatial discretizations mentioned above, the spatial error is $O(h^...


5

It could be that your optimization problem is very badly conditioned when using only the accelerometer data. In other words, the accelerometer data might not sufficiently constrain the parameters so that many different paths adequately fit the data. In terms of the minimization problem this means that you'd have a large "flat spot" at the minimum of ...


5

Strictly speaking, it doesn't make sense to speak of a Courant number for the Navier-Stokes equations at finite Reynolds number. The Courant number is defined in terms of a characteristic velocity, and solutions of parabolic equations (like Navier-Stokes) aren't described by characteristics. The CFL condition implies that no explicit solver with finite ...


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