Solving DAE with higher index

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))$

• You should add your equations so that we can see the problem. Commented Nov 4, 2015 at 0:20
• @fatemeh: As written, your equations are unreadable. Also, these equations are a system of PDEs. Please delete your comment and instead write the DAE system (that is to say, write out the Method of Lines discretization you are solving) in the body of your question post. Commented Nov 6, 2015 at 7:09
• @GeoMatt22 .thanks.i inserted equations in the question. Commented Nov 6, 2015 at 10:55

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.

• THANKS Geoff. Do you mean i should use one of these methods:1. Pantelides' algorithm 2. dummy derivative ? I dont understand what is your meaning about IDA. IS it a method to reduce index too? Can i use it now? regards Commented Nov 5, 2015 at 7:58
• IDA is a software package to solve index-1 DAEs or Hessenberg form index-2 DAEs. It is not an index reduction method. Unless you list the form of your DAE in your question post, I cannot say whether you can use IDA to solve your DAE. However, if you elect to reduce the index of your DAE, you should use either Pantelides' algorithm or the dummy derivative method. Commented Nov 6, 2015 at 7:08
• Mathematica's NDSolve already includes algorithms from Lawrance Livermore State Laboratory's Sundials package(computation.llnl.gov/casc/sundials/main.html). So IDA can be automatically used when the system is differential-algebraic system. Per my experience here(mathematica.stackexchange.com/questions/107463/…), the Sundials IDA only handles problems with DAE index up to 1. Commented Feb 17, 2016 at 12:16

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;
];