Решение дифференциального уравнения в python с odeint

Question

Я пытаюсь решить это дифференциальное уравнение R · (dq/dt) + (q/C) = Vi · sin (w · t), поэтому у меня есть этот код:

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 

Итак, у меня есть это, теперь я буду интегрировать его с odeint

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

Но есть что-то неправильно, я думаю. Моя главная цель - получить q (t), а затем разделить его на C, чтобы получить Vc.

Vc=p/C 
print Vc 

Answer 1

Есть несколько вопросов здесь:

Отдела в python2 как комментатор уже отмечался, не отбрасывать целые числа с плавающей точкой в ​​процессе деления. Таким образом, вы хотите, чтобы убедиться, что ваши параметры поплавки:

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

Теперь вам не нужно, за что вы делаете, чтобы эти параметры в качестве аргументов. Давайте начнем проще, и просто использовать производную их как глобальные переменные:

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

Далее, ваше пространство имеет только начальную и конечную точки, потому что п равно 1! Давайте установим ni выше, скажем, 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 

Вы заметите, что это все еще не очень хорошо. Ваш управляющий сигнал составляет 1000 Гц, поэтому вам понадобится еще больше очков или меньший диапазон данных. Давайте попробуем от 0 до 0,02, с 1000 точек, что дает нам 50 очков за колебания:

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 

Но это все-таки неловко нестабильны! Зачем? Поскольку большинство таких решателей имеют как относительные, так и абсолютные допуски ошибок. В numpy они устанавливаются по умолчанию примерно в 1,4e-8. Таким образом, абсолютный погрешность погрешности приближается к вашим нормальным значениям q/p. Вы хотите, чтобы либо использовать единицы, которые будут иметь ваши ценности будут больше (вероятно, лучшим решением), или установить абсолютный допуск к соответственно меньшей стоимости (простое решение, как показано здесь):

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)