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); 

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:

enter image description here

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;

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:

enter image description here

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

enter image description here


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

Is my ODE solver for nonlinear ODE:s correct?


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

  • $\begingroup$ So ODE45 is much better that my simple algorithm? I want have more control of the simulation such as delay etc. In ODE45, I cannot have delays. $\endgroup$ – Daniel Mårtensson May 27 '18 at 18:09

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.