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

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

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

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

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

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

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

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


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

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

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


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

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


5

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


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

A useful rule of thumb (though it is not always a sufficient condition for stability) is stability of the linearized scheme. Since you have a method of lines discretization, you can think of this geometrically as the condition that the eigenvalues of the jacobian of your spatial discretization, multiplied by the time step size, lie inside the region of ...


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

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 in space for instance), you get : $$ \frac{u_i^{n+1} - u_i^n}{\Delta t}+\frac{u_i^n-u_{i-1}^{n}}{\Delta x} = 0 ~~~~~~~~(2) $$ ...


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


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