# Solving a system of differential algebraic equations (DAE)

I am trying to solve a system of 8 differential algebraic equations, where equations 3 and 5 are differential equations and the rest are constraints which need to be satisfied. Also I only know the initial condition $y[0]=b/r$, so I am letting Mathematica find the consistent initial conditions.

Using NDSolve gives me an error:

NDSolveValue::index: The DAE solver failed at x = 0. The solver is intended for index 1 DAE systems and structural analysis indicates that the DAE index is 3. The option Method -> {"IndexReduction" -> Automatic} may be used to reduce the index of the system.

Since I have constraints (invariants) I also tried Method -> "Projection" but it gives me error:

NDSolve::initf: The initialization of the method NDSolveStateSpace failed.

How can I solve this problem using Projection methods or any other in Mathematica?

(*Defining parameters*)
param = {a -> 0.072, b -> 0.04, r -> 0.05, ρ -> 0.06, σ -> 0.1, τ ->  1/3, α -> 0.1};

(*Defining variables to be evaluated*)
vars = {y, RM, sxx, sy, mxx, my, R, λ};

(*Defining the system to be solved, differential equations and constraints (invariants)*)
in1 = RM[x] -(1-Exp[(my[x] + σ*sy[x] - 1/2*(σ + sy[x])^2)τ + (σ + sy[x])*InverseCDF[NormalDistribution[0, 1], α]*Sqrt[τ]]);
in2 = sxx[x] - (1/RM[x] - 1)*(σ + sy[x]);
eq3 = sy[x] == sxx[x]*D[y[x], x]/y[x];
in4 = mxx[x] -(  1/RM[x]*a/  y[x] - ρ + (1/RM[x] - 1)*(my[x] + σ*sy[x] - R[x]) + (1 - 1/RM[x])*(σ + sy[x])^2);
eq5 = my[x] == D[y[x], x]/y[x]*mxx[x] + 1/2*D[y[x], {x, 2}]/y[x]*(sxx[x])^2; in6 = (r*(1 - x) + ρ*x)*  y[x] -( a*x/RM[x] + (1 - x/RM[x])*b);
in7 = b/y[x] + my[x] + σ*sy[x] - R[x] - (1 - x/RM[x])/( 1 - x)*(σ + sy[x])^2;
in8 = a/y[x] + my[x] + σ*sy[x] - R[x] -( 1/RM[x]*(σ + sy[x])^2 + λ[x]*y[x]*RM[x]);

eqall = {in1==0, in2==0, eq3, in4==0, eq5, in6==0, in7==0, in8==0} /. param;

inv = {in1 == 0, in2 == 0, in4 == 0, in6 == 0, in7 == 0, in8 == 0} /. param;


Here is the solution part:

(*Finding consistent initial conditions*)

init = NDSolve[eqall, vars, {x, 0, 0}];
value = First[{y[x], y'[x], y''[x], RM[x], sxx[x], sy[x], mxx[x], my[x], R[x], λ[x]} /. init /. x -> 0]

(*Using NDSolve*)
sol = NDSolve[{eqall, y[0] == value[[1]], y'[0] == value[[2]], y''[0] == value[[3]], RM[0] == value[[4]], sxx[0] == value[[5]], sy[0] == value[[6]], mxx[0] == value[[7]], my[0] == value[[8]],  R[0] == value[[9]]}, vars, {x, 0, 1}];

(*Using NDSolve Projection Methods*)
Prosol = NDSolve[{eqall, start}, vars, {x, 0, 1}, Method -> {"Projection", Method -> "ExplicitRungeKutta",  "Invariants" -> inv}, InterpolationOrder -> Automatic]
`

• Welcome! I suggest the following: 1) As you receive help, try to give it too, by answering questions in your area of expertise. 2) Take the tour and check the faqs! 3) When you see good questions and answers, vote them up by clicking the gray triangles, because the credibility of the system is based on the reputation gained by users sharing their knowledge. Remember to accept the answer, if any, that solves your problem, by clicking the checkmark sign! – user9660 Oct 10 '16 at 9:17
• You did not input initial value for $Lambda$ in sol – Prasad Mani Oct 2 '17 at 6:13