# Solving an algebraic-differential equation

I cannot find the relevant part in the documentation. Maybe you have a quick solution to my problem.

## Background

Sometimes you want to propagate in time some system of ODEs, but not all the variables are needed. Think of a matrix Schrödinger equation. Let us say, a sum of squares of unknown functions is desired. One may still solve ODEs of interest, then build the desired observable. But this is also an ideal case for an algebraic-differential solver. One just need to propagate everything, but having only the observable as an output.

My problem is the following messages

NDSolve::mconly: For the method IDA, only machine real code is available. Unable to continue with complex values or beyond floating-point exceptions.
NDSolve::icfail: Unable to find initial conditions that satisfy the residual function within specified tolerances. Try giving initial conditions for both values and derivatives of the functions.


## Minimal working example

h[t_]:={{1,1,5},{1,2,I},{5,-I,3}}
o={{0,0,5},{0,0,0},{5,0,0}}
dim=3;
vars=y[#]&/@Range[dim];
varsT=y[#][t]&/@Range[dim];
y0={0,1,0};
r=NDSolve[Join[eqs,ics],vars,{t,0,10}]//First;


As the result we get some oscillating functions, an observable can be easily computed and plotted

ψ=vars/.r
oT=Sum[Conjugate[ψ[[i]][t]]o[[i,j]]ψ[[j]][t],{i,dim},{j,dim}];
Plot[oT,{t,0,10}]


I want to avoid post-processing to evaluate the observable oT.

## Minimal not working example

Instead of evaluating the observable after solving ODEs, let us propagate it along with the unknowns. We add one additional algebraic equation and the respective initial condition

eqA={τ[t]==Sum[Conjugate[y[i][t]]o[[i,j]]y[j][t],{i,dim},{j,dim}]};
icA={τ[0]==0};


However, this fails with the error messages above

NDSolve[Join[eqs,eqA,ics,icA],τ,{t,0,10}]


What can be done here? I am not even sure the initial condition icA for the auxiliary variable τ is needed.

• Related: mathematica.stackexchange.com/q/167368/1871 mathematica.stackexchange.com/q/184281/1871 In short, nobody has found an effective way to adjust the DAE solver in this site for so long AFAIK, and the best thing we can do seems to be avoiding DAE solver in these cases. – xzczd Aug 16 at 13:39
• Re NDSolve::mconly: One might try to model the complex arithmetic in terms of real numbers, but it might be easier to solve the ODEs and build the desired observable. Another approach is to differentiate the algebraic equation and solve it as an ODE. – Michael E2 Aug 16 at 13:44
• @xzczd Especially your 2nd link is relevant. I was not able to find it... So it is a general and yet unsolved problem. Should I delete the post to avoid duplication? – yarchik Aug 16 at 13:47
• I think it's OK to keep it, it contains a good minimal example. Perhaps someday NDSolve will be able to directly solve this, then we can come back and kill the unanswered question :) . – xzczd Aug 16 at 13:50