2015-03-19 49 views
0

我正在使用python ODE來求解Van der Pol方程。我想要做的是找到一次調用積分(t)方法期間調用ODE的派生數量的次數。我使用的解算器是isoda。我試圖找到在沒有jacbian函數的情況下在一次調用期間調用衍生物的次數,以及包含jacobian函數時調用jacobian函數的次數。lsoda.error:處理參數列表中的回調jac失敗

def modified_vanderpol_integrator_demo(dt_min=1e-3, dt_max=10., mu=100.): 
    # define a class to store the number of calls 
    class F: 

     def __init__(self): 
      self.callsydot = 0 
      self.callsjac = 0 

     def ydot(self,t, y, mu): 
      self.callsydot += 1 
      x, v = y[0], y[1] 
      return array([v, mu*(1- x*x)*v - x]) 
     def jac(self, t, y, mu): 
      self.callsjac += 1 

      return array([[0, 1], [2*mu*v*x - 1, mu*(1-x*x)]]) 

    Dt = linspace(dt_min, dt_max, 1000) 
    Num_of_times1 = [] 
    Num_of_times2 = [] 
    for dt in Dt: 
     xinitial = array([1.0, 2.0]) # initial value 
     f1 = F() 
     r1 = ode(f1.ydot) 
     r1.set_integrator('lsoda') 
     r1.set_initial_value(xinitial) 
     r1.set_f_params(mu)  
     r1.integrate(dt) 
     Num_of_times1.append(f1.callsydot) 

     f2 = F() 
     r2 = ode(f2.ydot, f2.jac) 
     r2.set_integrator('lsoda', with_jacobian=True) 
     r2.set_initial_value(xinitial) 
     r2.set_f_params(mu)  
     r2.integrate(dt) 
     Num_of_times2.append(f2.callsjac) 


    plt.plot(Dt, Num_of_times1) 
    plt.plot(Dt, Num_of_times2) 
    plt.show() 

當我運行該腳本,我得到的消息

create_cb_arglist: Failed to build argument list (siz) with enough arguments (tot-opt) required by user-supplied function (siz,tot,opt=2,3,0). 
Traceback (most recent call last): 
    File "stiffequations.py", line 118, in <module> 
    modified_vanderpol_integrator_demo(dt_min=1e-3, dt_max=10., mu=100.) 
    File "stiffequations.py", line 111, in modified_vanderpol_integrator_demo 
    r2.integrate(dt) 
    File "/usr/local/lib/python2.7/dist-packages/scipy/integrate/_ode.py", line 394, in integrate 
    self.f_params, self.jac_params) 
    File "/usr/local/lib/python2.7/dist-packages/scipy/integrate/_ode.py", line 1193, in run 
    y1, t, istate = self.runner(*args) 
lsoda.error: failed in processing argument list for call-back jac. 

爲什麼出現這種情況,如何解決? 謝謝。

回答

0

我發現我的腳本有什麼問題。我忘了在雅可比方程中設置參數。

r2.set_jac_params(mu)