I'm trying to solve a system of DAE (2 ODE and 1 algebraic equation) using the solver IDA from Sundials (https://computation.llnl.gov/projects/sundials/ida), through the Python package scikits.odes (https://scikits-odes.readthedocs.io).
I'm using scikits.odes 2.4.0, Sundials 3.1.1 and Python 3.6 64bit.
Here is the code :
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)
When I run it, I get the following error :
[IDA ERROR] IDASolve
At t = 0 and h = 2.86102e-07, the corrector convergence failed repeatedly or with |h| = hmin.
Do you have any idea of from where it could come from?
The immediate cause should be that the initial vector is not a consistent state, as it violates the algebraic part
0 == m1 @ y3 + v1
as y1=[0,0] and v1=[0.3, 9]*1e-4 is non-zero.
Apart from that, as far as I can see, your system has index 2, this requires a specialized DAE solver. DASPK, which is used in Sundials/IDA, generally only solves index-1 DAE. Depending on the version, certain special classes of index-2 DAE can also be solved. From the R wrapper documentation one can learn that if the maximal differentiation orders of the variables are known, that solutions for index up to 3 can be obtained. I do not know if this python wrapper is prepared for that.
The system
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
can be transformed into an ODE kernel and a state reconstruction procedure. The reduced state for the ODE consists of the components [ b, c3 ] with the equations
b' = -(c3 + v2)/m
c3' = 0.5*(f(c3)-f(v11+v12-c3)+v3-b)/M
and for the state reconstruction (starting from b,c3,c3' using the ODE function for derivatives)
c1 = v11+v12-c3
c2 = v22
a1 = f(c1)+b+M*c3'
a2 = f(c2)+a1
If one wanted to include the full state in the ODE one would need to differentiate all equations (except possibly the third) once more (thus differentiation index 2). Then the state of the DAE is obtained from the state of the ODE system (containing some derivative components) via projection.