CODEWITHSUNDEEP
odegekko
Chris Seeling
Chris Seeling

Reputation: 656

Simulation of Lorenz System crashes

I only just started exploring Gekko and tried simulating the Lorenz ODE system. Unfortunately, I get an error ("no solution found") for a simple std case the runs fine using scipy.

The problem solves fine if I only integrate up to say time=0.5 instead of 1.0

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

m = GEKKO()
m.time = np.arange(0.0, 1.0, 0.01)

sigma = 10.; rho   = 28.0; beta  = 8./3.

x = m.Var(value=10); y = m.Var(value=10); z = m.Var(value=10)
t = m.Param(value=m.time)

m.Equation(x.dt()== sigma*(y - x))
m.Equation(y.dt()== x*(rho -z) - y)
m.Equation(z.dt()== x*y - beta*z)
m.options.IMODE = 4
m.options.nodes = 4
m.solve(disp=False)

plt.plot(x.value, y.value)
plt.show()

Upvotes: 2

Views: 49

Answers (1)

John Hedengren
John Hedengren

Reputation: 14346

Sequential simulation m.options.IMODE=7 solves successfully. The simultaneous simulation is a larger problem (1782 Variables/Equations vs. 18 Variables/Equations). It also solves successfully if you reduce m.options.NODES=3 (1188 Variables) or increase the maximum iterations m.options.MAX_ITER=300. It failed previously because it needed 267 iterations to find the solution and the maximum iteration limit for IPOPT is 200 by default.

Lorenz

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

m = GEKKO()
m.time = np.arange(0.0, 1.0, 0.01)

sigma = 10.; rho   = 28.0; beta  = 8./3.

x = m.Var(value=10); y = m.Var(value=10); z = m.Var(value=10)
t = m.Param(value=m.time)

m.Equation(x.dt()== sigma*(y - x))
m.Equation(y.dt()== x*(rho -z) - y)
m.Equation(z.dt()== x*y - beta*z)
m.options.IMODE = 7
m.options.nodes = 4
m.solve(disp=True)

plt.plot(x.value, y.value)
plt.show()

Upvotes: 1

Related Questions