首页 文章

用odeint求解python中的微分方程[关闭]

提问于
浏览
-5

我试图解决微分方程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

1 回答

  • 2

    这里有一些问题:

    正如评论者所指出的那样,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 Hz,因此您需要更多的点或更小的数据范围 . 让我们尝试从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)
    

相关问题