Unstable simulation for Euler-ODE solver - but not with ODE45

Question

I have made an algorithm which can simulate both nonlinear ODE:s and linear ODE:s.

I begin with the linear.

Assume that we got a LTI-system:

$$\dot x(t) = Ax(t) + Bu(t) $$

And we turn it to discrete

$$x(k+1) = A_dx(k) + B_du(k)$$

Then we can simulate that by using this algorithm.

x = x0 % Initial state vector
for k = length(time)
  x = Ad*x + Bd*u(k); 
end

Because $x(k+1)$ is the next state vector $x(k)$

Assume that we have the system:

$$\dot x = \begin{bmatrix} 0 & 1 \\ -3 & -0.1 \end{bmatrix}\begin{bmatrix} x_1\\ x_2 \end{bmatrix} + \begin{bmatrix} 0\\ 1 \end{bmatrix}\begin{bmatrix} u_1 \end{bmatrix}$$

and the sample time $h = 0.10101$. The step answer simulation results look like this:

Where $y1$ is the output and $y2$ is the derivative of the output.

Now for the nonlinear one.

Assume that we have a NLTI-system described at this form:

$$\dot x = f(x, u)$$

and we separate the system at this form:

$$\dot x = A(x,u) + B(x,u)$$

Then we can simulate that system by using this algorithm:

for k = 1:size(t,2) 
  dxdt = A(x, u) + B(x, u);
  % update state
  x = x + h*dxdt;
end

Where $h$ is the sample time

Assume that we have the system:

$$\dot x = \begin{bmatrix} 0 & x_2 \\ -3x_1 & -0.1x_2^2 \end{bmatrix} + \begin{bmatrix} 0\\ u_1 \end{bmatrix}$$

and $h = 0.10101$. The step answer simulation results look like this:

But using ODE45 from Octave/MATLAB with the NLTI-system, the step answer simulation look like this:

Question:

Is runge-kutta method a much more stable ODE solver compared to Euler ODE solver?

Is my ODE solver for nonlinear ODE:s correct?

Answer 1

The Euler methods, interpreted as one-step methods, are Runge-Kutta methods.

The error order of the Euler methods is $1$, thus you get an error $O(t\cdot h)$ at time $t$, or more precisely for larger $t$ the error is proportional to $(e^{Lt}-1)h$ where $L$ is a Lipschitz constant of the ODE system. In your case, $L=3$ and with $t=10$ you get a possible error scale of $0.1\cdot e^{30}\approx 10^{12}$ which tells you that the Euler values are not very reliable on the considered interval.

ode45 uses the embedded order 4/5 Dormand-Prince method (not the classical RK4 method) with adaptive step size. This is as far away from the fixed-step Euler solver as is possible in Runge-Kutta solvers. You can expect that the default error tolerances ($10^{-6}$ or $10^{-9}$, check documentation) are actually true error estimates, so that on the level of the graphical representation the solution is exact..