Wrong solution for a simple ODE system

[ivlis]

I am trying to solve a simple ODE system:

du/dt = v
dv/dt = - u

with initial conditions u[0] = 1 and v[0] = 0 The anaytical solution is u[t] = cos(t) and v[t] = sin(t)

I am using the following GEKKO code:

LL = GEKKO()
LL.time = np.linspace(0,24*np.pi,1000)

S = [LL.Var(value=1), LL.Var(value=0)]

t = LL.Param(value=LL.time)

LL.Equations([S[0].dt() == -S[1], S[1].dt() == S[0]])

LL.options.IMODE = 4
LL.solve(disp=False)

The result is the following:

For some reason, the solution has an exponential decay....

Wolfram Mathematica solves the same equations just fine.

[APMonitor]

Some problems require more accuracy to avoid things like numerical drift. By default, GEKKO uses the lowest level of accuracy with NODES=2 (available options=2-6). Setting NODES>=3 as:

LL.options.NODES = 3

fixes the issue. Here is some documentation on NODES (order of polynomial used in Orthogonal Collocation on Finite Elements): https://apmonitor.com/wiki/index.php/Main/OptionApmNodes Here is the theory and examples on the method that GEKKO uses: http://apmonitor.com/do/index.php/Main/OrthogonalCollocation

from gekko import GEKKO
import numpy as np
import matplotlib.pyplot as plt

LL = GEKKO()
LL.time = np.linspace(0,24*np.pi,1000)

S = [LL.Var(value=1), LL.Var(value=0)]

t = LL.Param(value=LL.time)

LL.Equations([S[0].dt() == -S[1], S[1].dt() == S[0]])

# LL.Equation(sz.dt()== 0)

LL.options.IMODE = 4
LL.options.NODES = 3
LL.solve(disp=False)

plt.plot(LL.time,S[0].value,'b-')
plt.plot(LL.time,S[1].value,'r--')
plt.savefig('results.png')

We could switch the default to NODES=3 but then there are other models such as Autoregressive eXogenous Inputs (ARX) and discrete state space Linear Time Invariant models that require NODES=2 for mixed discrete/continuous models. Thoughts?

[ivlis]

Wonderful, thank you!

[APMonitor]

Thanks for highlighting this issue.