4

I want to solve a second order differential equation with GEKKO. In the documentation there is only an example that shows you how to solve a first order equation. I can't figure out how to write the second derivative of y to make it work.

This is the example from the documentation for a first order differential equation.

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

m = GEKKO()
m.time = np.linspace(0,20,100)
k = 10

y = m.Var(value=5.0)
t = m.Param(value=m.time)
m.Equation(k*y.dt()==-t*y)
m.options.IMODE = 4
m.solve(disp=False)

plt.plot(m.time,y.value)
plt.xlabel('time')
plt.ylabel('y')

plt.show()
John Hedengren
  • 12,068
  • 1
  • 21
  • 25
Samuel
  • 41
  • 3

1 Answers1

1

The first derivative can be declared as dy = y.dt() and the second derivative as ddy = dy.dt()

import numpy as np
import matplotlib.pyplot as plt


m = GEKKO()
m.time = np.linspace(0,20,100)
k = 10

y = m.Var(value = 5.0)
t = m.Param(value=m.time)
dy = m.Var(value = 0.0)
ddy = m.Var(value = -5/10)


m.Equations([
    k*y.dt()==-t*y,
    dy == y.dt(),
    ddy == dy.dt()

])

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


plt.figure()
plt.plot(m.time,y.value, label = 'y')
plt.xlabel('time')
plt.plot(m.time, dy.value, label  = 'dy')
plt.xlabel('time')
plt.plot(m.time, ddy.value, label  = 'ddy')
plt.xlabel('time')
plt.legend()
plt.show()

enter image description here

You can find more information here : https://apmonitor.com/wiki/index.php/Apps/2ndOrderDifferential

reyPanda
  • 510
  • 2
  • 11
  • 1
    To get a second derivative ddy, you'll want something like dy=y.dt(), ddy=dy.dt(). – aknight0 Jul 17 '19 at 21:49
  • well i get the concept but in your example the orange curce clearly is not the derivation of the blue one which should be the case, shouldnt it? So there must be smth wrong... And also a more direct way to solve this would be super nice! – Samuel Jul 18 '19 at 11:09
  • oh sorry i got confused by your notation it should be the other way around and then its right. – Samuel Jul 18 '19 at 11:14
  • oh ok. I totally misunderstood your question. I believe the edited answer is now correct. – reyPanda Jul 18 '19 at 15:24
  • The new solution is essentially what @aknight0 is proposing – reyPanda Jul 18 '19 at 15:30
  • Well what you suppose enabels me to get the second derivative of a function, which is nice. But i want to solve an equation that includes second order derivation: For example I start with tht problem: ddy = -k*t and from this i want to get y. Following @reyPanda I definde my equation like this `m.Equations([ dy == y.dt(), ddy == dy.dt(), ddy == -k*t ])` But when I try to solve it i get the error message `No JSON object could be decoded` and I have no idea what that means.... – Samuel Jul 22 '19 at 09:47