MATLABnumerical methodsordinary differential equations

Consider the following ODE system to be solved numerically:

\begin{align}

\dot{x}_1 &= -x_1

\\

\dot{x}_2 &= -2x_2 + \beta x_4

\\

\dot{x}_3 &= 2x_2

\\

\dot{x}_4 &= x_1 – \beta x_4

\end{align}

When I solved the system numerically, I found that increasing the value of $\beta$ increases the number of iterations until convergence in some solvers (e.g. Matlab ode113), but the number of iterations stays relatively constant in other solvers (e.g. Matlab ode15s), and in fact it does it in fewer iterations.

My question: how does the magnitude of $\beta$ in this system affect the rate of convergence of a numerical method? Is there any specific relation that clearly shows that?

### ode15s

### ode113

#### Best Answer

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 solver ode113** 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.

#### Related Solutions

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

###### Related Question

- Implicit method to solve a system of ODEs
- Mass matrix second order differential equation to first order ODEs
- Finite difference method: numerical solution is not approximate the exact solution.
- Implicit numerical ODE solvers: still unconditionally stable if Newton iterations limited