22

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

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


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

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


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


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

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

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

You solve the 1D-advection equation with $c$ a constant velocity: $$\frac{\partial u}{\partial t} + c\frac{\partial u}{\partial x}=0~~~~~~~~(1)$$ When you discretize this equation (with an explicit scheme in time and an upwind scheme for a positive velocity $c \geq 0$ in space for instance), you get: $$\frac{u_i^{n+1} - u_i^n}{\Delta t}+c\frac{u_i^n-u_{i-1}^{...


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

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


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

You can rewrite your finite difference method into the form \begin{align}\left[\begin{matrix} T^{t+1} \\ h^{t+1}\end{matrix} \right]=\left[\begin{matrix}A_{11} & A_{12} \\ A_{21} & A_{22} \end{matrix} \right] \left[\begin{matrix} T^n\\h^n\end{matrix} \right] = \left[\begin{matrix}A_{11} & A_{12} \\ A_{21} & A_{22} \end{matrix} \right]^n\...


6

You can define a method/algorithm stable if during the various steps it not amplifies excessively the errors on the data. So you introduce some bounded estimate. Note that this kind of errors, monitored by the stability check, are different from an error occurred by the ill-condition of a problem, i.e. in the stability you check the errors propagated by ...


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

(My apologies for the terrible typesetting) Let me start with the definition of backward stability of an algorithm: An algorithm $A$ to solve a problem $P$ is called backward stable, if for all $x$ in the input space there exists $y$ such that $||y-x||/||x||\approx \epsilon_m$ and $A(x) = P(y)$. Trefethen & Bau summarizes this as "Backward stable ...


5

In the following, I am basically rephrasing p. 24 of "MCMC using Hamiltonian Dynamics" by Radford Neal. At least for leapfrog integration, one can analytically calculate the maximum time step for quadratic Hamiltonians. For such systems, a leapfrog step is a linear mapping $(q(t), p(t) \mapsto (q(t+\epsilon), p(t+\epsilon))$. Stability then depends on the ...


5

In the advection dominated case the problem can develop physics which is invisible to a computational mesh that is too coarse (say by having elements which contain many wavelengths). This leads to instability because neighboring elements could disagree considerably on what they believe are the right values to represent this physics (think interpolating high ...


5

Analytical problem : what you are expecting is positive diffusion : you want the $T_i$ values to spread over your domain as time passes to eventually reach $T_i(t\rightarrow \infty) = cte$ if $\alpha$ was a negative number, you would have what is called negative diffusion : you'll have exactly the opposite i.e. the gradients will get greater through time. ...


5

There are two different flavors of smoothing and stability. The spurious oscillations in convection-diffusion problems are not an artifact of the linear solver but an inevitable artifact of the discretization method. Any linear solver for non-symmetric matrices you use will give you the same wiggly answer, or it's wrong. Multi-grid is one such linear solver. ...


5

Eigenvalues with zero eigenvalue correspond to purely oscillatory modes. You can see it by diagonalising the system. Your matrix $A$ can be written as \begin{equation} A = P \Sigma P^{-1} \end{equation} where $\Sigma$ is a diagonal matrix with entries $\lambda_n$ corresponding to the eigenvalues of $Q$. You can now transform your ODE to \begin{equation} dQ/...


5

This is comment, but too long for comment so I write here. I think is better you edit the post with more information, now I explain why. Here you can find some general advice to improve precision. The problem you find I think it's possible that it come from the floating point representation and the relative effect in the machine operation. When you do an ...


5

There is something very basic that you should know about hyperbolic problems. Consider the most basic example $\partial_tu+a\partial_xu=0$ with a numerical marching scheme of the form $$u_j^{n+1}=\sum_kc_ku_{j+k}^n$$. This covers all explicit schemes, and all implicit schemes like Crank Nicolson also, if you begin by solving the tridiagonal system. It ...


5

Each Lagrange multiplier corresponds to a constraint. So if the space of Lagrange multipliers is too large, then you have too many constraints that can no longer be all satisfied at the same time without significantly restricting the number of unknowns you have available to satisfy the physics of the problem. This is what happens in locking: each constraint ...


5

The saddle point matrix you must solve takes the following form: $$\begin{bmatrix} A & B^T \\ B \end{bmatrix},$$ where $A$ is the unconstrained matrix and $B$ is the matrix that observes the values of the solution at nodes along the immersed interface. The conditioning of the saddle point matrix is, in part, controlled by the minimum singular value of $...


5

Orthogonal matrices are about as well-conditioned as you can get, but numerical errors still occur. One common error is loss of orthogonality. A fix for this could be to re-orthogonalize your columns after some number of multiplications. You can do this by just taking the QR decomposition of your matrix after some number of products and taking the orthogonal ...


5

So, you want to invert your matrix $A=\Phi^T\Phi$. For $A$ to be invertible it must not have zero eigenvalues. We can show that $A$ is positive semi-definite as follows. Positive semi-definite means that the eigenvalues of $A$ are $\geq 0$. This is equivalent to showing $y^TAy \geq 0, \forall y \neq 0 $. $$ y^TAy = y^T\Phi^T\Phi{y}=(\Phi{y})^T(\Phi{y}) \geq ...


5

As hardmath mentioned in the comments above, if data space $D$ is empty, then the first half of the condition "for all $X\in D$ there exist $Y\in D$ such that $||X-Y||/||X||\approx\epsilon_m$ and $A(X)=P(Y)$" becomes vacuously true. So for backward stability analysis you only need to show the second part. In this case, it is simply proving that the ...


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