Как исправить "IDASolve: сбой корректора" при решении DAE с scikits.odes и Sundial?

Я пытаюсь решить систему DAE (2 ODE и 1 алгебраическое уравнение), используя решатель IDA из Солнечных часов ( https://computation.llnl.gov/projects/sundials/ida), через пакет Python scikits.odes ( https://scikits-odes.readthedocs.io/).

Я использую scikits.odes 2.4.0, Солнечные часы 3.1.1 и Python 3.6 64bit.

Вот код:

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)

Когда я запускаю его, я получаю следующую ошибку:

[IDA ERROR]  IDASolve
  At t = 0 and h = 2.86102e-07, the corrector convergence failed repeatedly or with |h| = hmin.

Есть ли у вас какие-либо идеи о том, откуда это может прийти?

Источник

2 ответа

Решение

Непосредственной причиной должно быть то, что начальный вектор не является непротиворечивым состоянием, так как он нарушает алгебраическую часть

0 == m1 @ y3 + v1

как y1=[0,0] а также v1=[0.3, 9]*1e-4 ненулевой

Кроме того, насколько я вижу, ваша система имеет индекс 2, для этого требуется специализированный решатель DAE. DASPK, который используется в Sundials/IDA, обычно решает только DAE индекса-1. В зависимости от версии, определенные специальные классы DAE индекса-2 также могут быть решены. Из документации оболочки R можно узнать, что, если известны максимальные порядки дифференцирования переменных, можно получить решения для индекса до 3. Я не знаю, подготовлена ​​ли эта оболочка Python для этого.

Решение без решателя DAE с помощью ручного уменьшения индекса

Система

0 = -c1+c2-c3 + v11
0 =    -c2    + v12

m*b' = -c1 - v2

M*c1' = f(c1) - a1     + b 
M*c2' = f(c2) + a1-a2   
M*c3' = f(c3) - a1     + v3

может быть преобразовано в ядро ​​ODE и процедуру восстановления состояния. Сокращенное состояние для ODE состоит из компонентов [ b, c3 ] с уравнениями

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

и для государственной реконструкции (начиная с b,c3,c3' используя функцию ODE для производных)

c1 = v11+v12-c3
c2 = v22
a1 = f(c1)+b+M*c3'
a2 = f(c2)+a1

Если кто-то хочет включить полное состояние в ODE, ему нужно будет еще раз дифференцировать все уравнения (кроме, возможно, третьего) (таким образом, индекс дифференциации 2). Затем состояние DAE получается из состояния системы ODE (содержащей некоторые производные компоненты) посредством проекции.

Источник

LutzL прав насчет начальных условий. Спасибо также за ответ ACHindmarsh на форуме SUNDIALS ( http://sundials.2283335.n4.nabble.com/IDA-ERROR-IDASolve-the-corrector-convergence-failed-repeatedly-or-with-h-hmin-td4655649.html), я более подробно изучил документацию Scikits.Odes ( https://github.com/bmcage/odes/blob/master/scikits/odes/sundials/ida.pyx) и обнаружил, что Оболочка для решателя IDA может принять вариант compute_initcond to указывает начальное условие, которое мы хотим, чтобы решатель вычислял сам. Таким образом, я установил эту опцию на 'y0' и решатель теперь успешно интегрирует мою систему.

Источник
Другие вопросы по тегам