Solving DAE with higher index

Question

I want to solve 7 coupled differential equations in MATLAB. I used the method of lines (MOL) to convert them to a system of DAEs. I then attempted to solve this system using ode15s. Unfortunately an error occurred that said: the index of equation is greater than 1. I searched and know I should reduce the index of equation. But I dont know how can I do this in MATLAB?
I have simplified equations and assume T=CONSTANT so these  are 3 equations:

$1:'(∂C_(g))/∂t+ u_g  (∂C_(g))/∂z= -k_g A_o (C_(g)-C_(p) |_(r=r_OD ) )+A_o ε_f u_cp |_(r_OD ) C_(p) |_(r_OD )'$
$2:(∂C_(p))/∂t+(u_cp C_(p))/r+u_cp  (∂C_(p))/∂r-D_eff ((∂^2 C_(p))/(∂r^2 )+1/r  (∂C_(p))/∂r)=-((1-ε_f))/ε_f  ρ_f  (∂q_)/∂t$
$3:(∂q)/∂t=ks (q^eq-q)$
$4:q^eq=(qs.b.P_)/〖(1+(b.P_ )^n)〗^(1/n)$
$5:P=C(p)RT$

and these are initial and boundry conditions:
at $t=0$  $q=0$, $cp=0$ and $cg=0$
at $z=0$ $cg=constant$
ar $r=riD$  $(∂C_(i,p))/∂r=0$
at $r=roD$   $Deff.(∂C_(i,p))/∂r=kg.(cg(z)-cp(r=roD)$

I convert them to system of odes in time.
these are equations:
$i$ is used for lenght(z).and $j$ is used for radius(r)

$1':dcg(i,j)=((-1).*(Ug).*(((cg(i,j))-(cg(i-1,j)))./((deltaz))))-((kg)*(Ao).*((cg(i,nr2))-(cp(i,nr2))))+((Ao)*(epsilonf).*(Ucp)*(cp(i,nr2)))$
$2':dcp(i,j)=(((-1).*(Ucp).*(cp(i,j)))./((Ri)+((j-1)*(deltar))))-((Ucp).*(((cp(i,j))-(cp(i,j-1)))./((deltar))))+((Deff).*((((cp(i,j+1)))-(2*cp(i,j))+(cp(i,j-1)))./(deltar^2))+(((1)/((Ri)+((j-1)*(deltar)))).*(((cp(i,j))-(cp(i,j-1)))./((deltar)))))-((((1)-(epsilonf))/(epsilonf))*(ROf).*(dq(i,j)))$
$3': dq(i,j)=ks*(qeq(i,j)-q(i,j))$

Geoff Oxberry · Accepted Answer

There are many definitions of DAE index. MATLAB is probably referring to the differentiation index. A few integrators can solve Hessenberg form index-2 DAEs (I believe IDA can do this, along with a few other packages), but most require the index to be 1 or less. Reducing the index of your DAE requires manipulating the underlying equations in the DAE, e.g., using Pantelides' algorithm or the dummy derivative method. Usually, this sort of thing is done in software that has some sort of algebraic manipulation capabilities, like Mathematica, Maple, Sage, or SymPy, so in your case, you would use one of the above approaches to derive an equivalent reduced-index DAE, and then transcribe it into MATLAB and solve it with ode15s.

Correct answer by Geoff Oxberry on May 9, 2021

Chris Rackauckas · Answer

If you use Julia's DifferentialEquations.jl it can automatically fix the index of your equations. A tutorial of this is shown in ModelingToolkit.jl. For example, we can write down an index-3 DAE:
using ModelingToolkit
using LinearAlgebra
using OrdinaryDiffEq

function pendulum!(du, u, p, t)
    x, dx, y, dy, T = u
    g, L = p
    du[1] = dx
    du[2] = T*x
    du[3] = dy
    du[4] = T*y - g
    du[5] = x^2 + y^2 - L^2
    return nothing
end
pendulum_fun! = ODEFunction(pendulum!, mass_matrix=Diagonal([1,1,1,1,0]))
u0 = [1.0, 0, 0, 0, 0]
p = [9.8, 1]
tspan = (0, 10.0)
pendulum_prob = ODEProblem(pendulum_fun!, u0, tspan, p)

Then we can use modelingtoolkitize to get the symbolic form of the equations, and then dae_index_lowering to receive the index-1 form:
traced_sys = modelingtoolkitize(pendulum_prob)
pendulum_sys = structural_simplify(dae_index_lowering(traced_sys))

and now you can solve this system with an index-1 DAE solver:
prob = ODAEProblem(pendulum_sys, Pair[], tspan)
sol = solve(prob, Tsit5(),abstol=1e-8,reltol=1e-8)

using Plots
plot(sol, vars=states(traced_sys))

Now, if you find the need to decelerate your code, you can then use build_function directly to transform the index-1 result into MATLAB code.
build_function(map(x->x.rhs,equations(pendulum_sys)),
                [independent_variable(pendulum_sys)],
                states(pendulum_sys),
                parameters(pendulum_sys),
                target=ModelingToolkit.MATLABTarget()
                )

which results in:
diffeqf = @(t,internal_var___u) [
  internal_var___u(2);
  internal_var___u(1) * internal_var___u(5);
  internal_var___u(4);
  -1 * internal_var___p(1) + internal_var___u(3) * internal_var___u(5);
  2 * internal_var___u(2) ^ 2 + 2 * internal_var___u(4) ^ 2 + 2 * internal_var___u(3) * (-1 * internal_var___p(1) + internal_var___u(3) * internal_var___u(5)) + 2 * internal_var___u(5) * internal_var___u(1) ^ 2;
];

Solving DAE with higher index

2 Answers

Add your own answers!

Ask a Question