# Tag Info

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

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 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 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 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 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 Bluntly speaking, SUPG and alike and RANS are different approaches to different problems that, however, have the same name - instability - and the same phenomenology - the failure of numerical routines. RANS is used to cope with turbulence as an instability of the equation. If a flow is or becomes turbulent the describing equations are instable, e.g. ... 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 I have little to contribute here other than to point out that whenever numerical methods have trouble with hyperbolic equations (and converge to the wrong solution), it isn't usually because of shocks. Rather, the areas they have difficulty with it are rarefaction waves -- where the solution is smooth. Another example of something that appears to be ... 5 Look at the cell Péclet number$$\mathrm{Pe}_h = h v / K$$where$h$is mesh size,$v$is the magnitude of velocity, and$K$is diffusivity. It is analogous to cell Reynolds number for the momentum equation and is small when "thermal diffusivity is large compared to advection". It is common common in macro-scale fluid dynamics that thermal diffusivity$K\$...

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