How does a constant magnitude in a system of ODEs affect the convergence of a numerical method – Math Solves Everything (2024)

Your system converges to a fixed point with $x_1^*=x_3^*=x_4^*=0$ and $x_2^*=c_1+c_2+c_3+c_4$ where $x(0)=c$. That last follows from the fact that the sum over the derivatives cancels all the terms on the right side, so that the sum over the components is constant.

The non-stiff solverode113 uses Adams-Bashford-Moulton, but in a PECE implementation which makes this an explicit method. As such it has a bounded stability region for any order that is used internally. Assume that the order stabilizes. Then the step size controller will increase the step size $h$ until that boundary is reached and surpassed for $\lambda h$ where $\lambda=-\beta$ is the dominant eigenvalue of the Jacobian. Remember that the stability region is where the step propagator is contracting (non-expanding) if the flow around the exact solution is contracting (non-expanding).

At this point, outside the stability region, the distance from the fixed point will grow, and with it the step error, so that the step size gets regulated down. Good step size controllers will minimize the resulting oscillation. The end result still is that for the most part of the integration interval you get $\beta h=const.$ where the constant is some value typically in the interval $[2,3]$. The number of evaluations over some interval $[0,T]$ will then be about $N=T/h$ which is proportional to $\beta$.

Note that this behavior is largely independent of the target error tolerance, the equilibrium step size will not change, only the equilibrium distance from the stationary point will change.

In the stiff solver ode15s the BDF methods that are truly implemented as implicit methods do not have (all) such restricted stability regions. Thus the error estimator is less wrong and keeps the step size in the region that is optimal for the error tolerance. As the components that have reached their asymptotic value fall below the absolute error tolerance, they do not influence the step size controller. Thus the step size dynamic is largely oriented to approximate the $e^{-t}$ and $e^{-2t}$ components correctly, which is independent of $\beta$ and gives thus the almost constant step counts.

[Math] Rate of Convergence for Trapezoidal Method-System of Linear ODEs

Your exact value is in the second component not exact enough. With N ranging from 32 to 1024, one would expect for an order 2 method and a tame ODE system errors from 1e-3 to 1e-6. However, if the supposedly exact value has an error of magnitude 1e-4 to 1e-5, then the integration error at the last iterations get unnecessarily perturbed.

\begin{align} u_1(t) &= \exp(t)\\ (u_2(t)·\exp(t))' &= u_1(t)·\exp(t) = \exp(2·t)\\ u_2(t)·\exp(t) &= C + 0.5·\exp(2·t) \implies C = 0.5\\ u_2(t) &= (\exp(t)+\exp(-t))/2 = \cosh(t)\end{align}

Numerical Methods – Numerical Solution to System of ODEs Coupled to a PDE

The usual trick for handling such coupled PDE-ODE systems with different time scales is splitting. If we consider the vector of ten unknowns $\mathbf{U} = (E_p,p_{11}\dots p_{33})^\top$, then the system of differential equations can be written$$\frac{\partial}{\partial t} \mathbf{U} + \mathbf{A} \frac{\partial}{\partial z} \mathbf{U} = \mathbf{R}(\mathbf{U}) \, ,$$where $\mathbf{A}$ is a $10\times 10$ matrix with only nonzero coefficient $\mathbf{A}_{11} = c$,and\begin{aligned}\mathbf{R}(\mathbf{U}) &= \left(i k c p_{13},\dot{p}_{11}\dots \dot{p}_{33}\right)^\top\\&= \left(i k c p_{13},E_p(p_{13} - p_{31})\dots E_p(p_{13} + p_{31}) + E_c (p_{32}-p_{23})\right)^\top .\end{aligned}In the present case, the two involved physical time scales relate to the wave speed $c$, and to the relaxation spectrum deduced from $\mathbf{R}(\mathbf{U})$. Splitting consists in solving successively$$\frac{\partial}{\partial t} \mathbf{U} + \mathbf{A} \frac{\partial}{\partial z} \mathbf{U} = \mathbf{0}\, , \qquad (\mathbf{H}_a)$$and$$\frac{\partial}{\partial t} \mathbf{U} = \mathbf{R}(\mathbf{U})\, , \qquad\qquad (\mathbf{H}_b) $$with a dedicated method for each operator.For instance, the Strang splitting scheme writes$$\mathbf{U}_i^{n+1} = \mathbf{H}_b(\Delta t/2)\circ \mathbf{H}_a(\Delta t) \circ \mathbf{H}_b(\Delta t/2)\, \mathbf{U}_i^n \, ,$$where $\mathbf{U}_i^n$ approximates $\mathbf{U}(z_i,t_n)$ and $\Delta t$ is the time step.This scheme is second-order accurate if both suboperators are stable and second-order accurate. Possible dedicated methods for each operator are

• a second-order accurate solver for the linear transport equations $\mathbf{H}_a$, e.g. a Lax-Wendroff or ADER scheme. Such a scheme is stable under the CFL condition.

• a second-order accurate solver for the nonlinear ODE $\mathbf{H}_b$, e.g. a Runge-Kutta method. If $\mathbf{H}_b$ is stiff, then the stability condition for such a solver can be penalizing.