如何修复在使用化学试剂盒日晷呢？

Question

我试图通过Python包使用sundalials（https://computation.llnl.gov/projects/sundials/ida）中的求解器IDA来求解DAE系统（2个ODE和1个代数方程）化学试剂盒（https://scikits-odes.readthedocs.io）。在

我在用化学试剂盒2.4.0、日晷3.1.1和Python 3.6 64位。在

代码如下：

import numpy as np
from scikits.odes.dae import dae

SOLVER = 'ida'
extra_options = {'old_api': False, 'algebraic_vars_idx': [0, 1]}

##### Parameters #####

# vectors
v1 = np.array([3.e-05, 9.e-04])
v2 = np.array([-0.01])
v3 = np.array([100])

# matrices
m1 = np.array([[-1, 1, -1], [0, -1, 0]])
m2 = np.array([[1, 0, 0]])
m3 = np.array([[0, 0, 1]])
m4 = np.array([[0., 0., 0., 0., 0., 0.],
               [0., 0., 0., 0., 0., 0.],
               [0., 0., 2000., 0., 0., 0.],
               [0., 0., 0., 13e3, 0., 0.],
               [0., 0., 0., 0., 13e3, 0.],
               [0., 0., 0., 0., 0., 13e3]])

##### Equations #####

def f(_, y):
    y1 = y[:2]
    y2 = y[2:3]
    y3 = y[3:]
    return np.hstack((m1 @ y3 + v1,
                      - m2 @ y3 - v2,
                      - 2e2 * np.abs(y3*1000) ** 0.852 * y3 + m1.T @ y1 + m2.T @ y2 + m3.T @ v3))

def right_hand_side(_, y, ydot, residue):

    residue[:] = f(_, y) - m4 @ ydot

    return 0

##### Initial conditions and time grid #####

tout = np.array([0.,  300.])

y_initial = np.array([0., 0., 0., 0., 0., 0.])

ydot_initial = np.array([0., 0., 0., 0., 0., 0.])

##### Solver #####

dae_solver = dae(SOLVER, right_hand_side, **extra_options)
output = dae_solver.solve(tout, y_initial, ydot_initial)
print(output.values.y)

当我运行它时，我得到以下错误：

^{pr2}$

你知道它是从哪里来的吗？在

Answer 1

LutzL对初始条件的看法是正确的。多亏了A.C.Hindmarsh在日晷论坛（http://sundials.2283335.n4.nabble.com/IDA-ERROR-IDASolve-the-corrector-convergence-failed-repeatedly-or-with-h-hmin-td4655649.html）上的回答，我对化学试剂盒（https://github.com/bmcage/odes/blob/master/scikits/odes/sundials/ida.pyx），我发现求解器IDA的包装器可以使用一个选项compute_initcond来指定我们希望解算器自己计算的初始条件。因此，我将这个选项设置为'y0'，解算器现在成功地集成了我的系统。在

Answer 2

直接原因应该是初始向量不是一致状态，因为它违反了代数部分

0 == m1 @ y3 + v1

因为y1=[0,0]和{}是非零的。在

除此之外，据我所知，您的系统有索引2，这需要一个专门的DAE解算器。DASPK用于日晷/IDA，通常只解决index-1dae。根据版本的不同，还可以解决index-2 DAE的某些特殊类。从R包装器文档中我们可以了解到，如果变量的最大微分阶数已知，则可以得到索引高达3的解。我不知道这个python包装器是否准备好了。在

---

无需DAE解算器的解决方案通过手动索引缩减

系统

^{pr2}$

可以转换成ODE内核和状态重建过程。ODE的约化状态由[ b, c3 ]和方程组成

b'  = -(c3 + v2)/m
c3' = 0.5*(f(c3)-f(v11+v12-c3)+v3-b)/M

对于状态重构（从b,c3,c3'开始，使用微分的ODE函数）

^{4}$

如果要在ODE中包含完整的状态，则需要再次微分所有方程（可能除了第三个方程）（因此微分指数2）。然后通过投影从ODE系统的状态（包含一些派生组件）中获得DAE的状态。在