0

I'm working with an ODE of the form:

a*dv/dt + (b+k1)*v + c*integral_0->t_(vdt) = k1*v1 + k2*integral_0->t_(v1dt)

I'm trying to implement odeint to get a solution for this system, but I'm not sure how to do that with an integral in the ODE. v1 is a known input, so that integral on the right side isn't a concern.

Spuds
  • 148
  • 1
  • 2
  • 13
  • I think you could derive this equation again for t, and solve the second order ode – xdze2 Aug 30 '18 at 08:12
  • 1
    Is this a correct interpretation of your equation? http://quicklatex.com/cache3/ae/ql_b8a0fa23c1efd633ffc409a0a3022aae_l3.png If so, you could consider including this image in your post instead. – tBuLi Aug 30 '18 at 09:26

1 Answers1

1

Set x as the integral of v, so that x'=v, x''=v', similarly x1 for v1, so that your equation reads as second order differential equation

a*x''+(b+k1)*x'+c*x=k1*v1+k2*x1

Given v1 as input, the state vector needs to contain the three integrated variables x, v, x1 which gives the ODE function

def odesys(y,t):
    x, v, x1 = y
    v1 = eval_v1(t)
    return [ v, (k1*v1+k2*x1 - (b+k1)*v-c*x )/a, v1 ]

To use it with odeint you can for example do

t = np.linspace(0,T,2001); # define the end time T before
y0 = [ 0, 0, 0 ]           # standard convention is that everything is zero for negative times
y = odeint(odesys, y0, t)  # add arguments for higher accuracy if needed
x, v, x1 = y.T             # transpose of a list of tuples is a tuple of lists
plt.plot(t,x); plt.show()  # as example that should work
Lutz Lehmann
  • 25,219
  • 2
  • 22
  • 51
  • So what would the odeint call look like? I'm not sure how to set up the y input for it. How do I define x and v? Sorry, this is my first time using odeint and I'm not finding the documentation very clear. – Spuds Aug 30 '18 at 18:57