-6

I am trying to solve that differential equation R·(dq/dt)+(q/C)=Vi·sin(w·t), so i have the this code:

import numpy as np
from numpy import *
import matplotlib.pyplot as plt
from math import pi, sin 
from scipy.integrate import odeint
C=10e-9
R=1000 #Ohmios
Vi=10 #V
w=2*pi*1000 #Hz
fc=1/(2*pi*R*C)
print fc
def der (q,t,args=(C,R,Vi,w)):
     I=((Vi/R)*sin(w*t))-(q/(R*C))
     return I

Ok, so i have this, now i am going to integrate it with odeint

a, b, ni = 0.0, 10.0, 1
tp=np.linspace(a, b, ni+1)
p=odeint(der, 0.0, tp)
print p

But there is something wrong i think. My main goal is to get q(t) and then divide it by C to get Vc.

Vc=p/C
print Vc
DiegoAlvrz
  • 1
  • 4
  • 3
    Please explain what is wrong. If there is an error, post it. – adarsh Apr 10 '15 at 20:49
  • This is python 2, right? R and Vi are integers, so Vi/R is an integer division so is 0, and so I is always zero. Make R=1000., Vi=10., and increase your ni. I also feel confident that your default argument for args isn't doing what you want, but that's a different issue. – Jonathan Dursi Apr 10 '15 at 20:51
  • When i begin with odeint, i don't know if i am doing it right, maybe i should add some arguments here? p=odeint(der, 0.0, tp, args=(w, R, C, Vi)? but it says ''der() takes at most 3 arguments (6 given)'' – DiegoAlvrz Apr 10 '15 at 20:54
  • Yes it is @Jonathan Dursi i tryed that but i think it still wrong – DiegoAlvrz Apr 10 '15 at 20:59

1 Answers1

2

There are a few issues here:

Division in python2, as commenters have noted, doesn't cast integers to floats during division. So you'll want to make sure your parameters are floats:

 C=10e-9
 R=1000.0 #Ohmios
 Vi=10.0 #V
 w=2*pi*1000.0 #Hz

Now, you don't really need, for what you're doing, to have those parameters as arguments. Let's start simpler, and just have the derivative use these as global variables:

def der (q,t):
    I=((Vi/R)*sin(w*t))-(q/(R*C))
    return I

Next, your space only has the start and end points, because ni is 1! Let's set ni higher, to say, 1000:

a, b, ni = 0.0, 10.0, 1000
tp=np.linspace(a, b, ni+1)
p=odeint(der, 0, tp)
plt.plot(tp,p) # let's plot instead of print; print p

You'll notice this still doesn't work very well. Your driving signal is 1000 Hz, so you'll need even more points, or a smaller data range. Let's try from 0 to 0.02, with 1000 points, giving us 50 points per oscillation:

a, b, ni = 0.0, 0.02, 1000
tp=np.linspace(a, b, ni+1)
p=odeint(der, 0, tp)
plt.plot(tp,p) # let's plot instead of print; print p

But this is still awkwardly unstable! Why? Because most ode solvers like this have both relative and absolute error tolerances. In numpy, these are both set by default to around 1.4e-8. Thus the absolute error tolerance is frighteningly close to your normal q/p values. You'll want to either use units that will have your values be larger (probably the better solution), or set the absolute tolerance to a correspondingly smaller value (the easier solution, shown here):

a, b, ni = 0.0, 0.02, 1000
tp=np.linspace(a, b, ni+1)
p=odeint(der, 0, tp, atol=1e-16)
plt.plot(tp,p)
cge
  • 9,552
  • 3
  • 32
  • 51