What makes a system be differential-algebraic equations?

Question

I am trying to solve the equations in the following way:

ndsovletmin = 0.;
ndsolvetmax = 100.;
EoMphi := {phi''[t] + 3. H[t] phi'[t] + phi[t] == 0.};
EoMa := {a'[t] == a[t] H[t], 3. H[t]^2 == phi'[t]^2/2. + phi[t]^2/2., Ne'[t] == H[t]};
bbc := {a[ndsovletmin] == 1., Ne[ndsovletmin] == 0., 
        phi[ndsovletmin] == 2., phi'[ndsovletmin] == 0.};
sol = NDSolve[{EoMphi, EoMa, bbc}, {phi[t], phi'[t], a[t], Ne[t]}, {t,
  ndsovletmin, ndsolvetmax}];

It gets a warning "NDSolve::pdord: Some of the functions have zero differential order, so the equations will be solved as a system of differential-algebraic equations.". But if I do a redefinition dphi[t] == phi'[t], ddphi[t] == dphi'[t], the warning disappears:

ndsovletmin = 0.;
ndsolvetmax = 100.;
EoMphi := {ddphi[t] + 3. H[t] dphi[t] + phi[t] == 0., 
           dphi[t] == phi'[t], ddphi[t] == dphi'[t]};
EoMa := {a'[t] == a[t] H[t], 3. H[t]^2 == dphi[t]^2/2. + phi[t]^2/2., Ne'[t] == H[t]};
bbc := {a[ndsovletmin] == 1., Ne[ndsovletmin] == 0., 
        phi[ndsovletmin] == 2., dphi[ndsovletmin] == 0.};
sol = NDSolve[{EoMphi, EoMa, bbc}, {phi[t], phi'[t], a[t], Ne[t]}, {t,
  ndsovletmin, ndsolvetmax}];

What's the difference?

You don't solve for H[t] from H'[t]. So H[t] has zero differential order. It seems that H[t] is positive. So try defining H[t_] := Sqrt[(phi'[t]^2 + phi[t]^2)/6.]; outside your system and then just use H[t]. When I try that I get no errors or warnings with your original method. Check all this very carefully to try to make certain this is correct. — Bill, Dec 16 '20 at 05:29
@Bill, I try what you say, but the warning still exist, I have edit my question and post the figure there. — JieJiang, Dec 17 '20 at 01:09

xzczd · Answer 1 · 2020-12-19T05:23:31.447

Both of the systems are DAE system, and NDSolve has used DAE solver to solve both of them. The difference is, in 2nd case, the system involves explicit algebraic equations 3. H[t]^2 == dphi[t]^2/2. + phi[t]^2/2. and ddphi[t] + 3. dphi[t] H[t] + phi[t] == 0. so NDSolve detects it's a DAE system immediately, while in 1st case, though the system actually doesn't involve derivative of H[t], it's not immediately obvious in NDSolve's eyes, so it first trys to transform the system to the standard form of ODE system (see this post for more discussion) and fails, thus spits out the pdord warning and turns to the DAE solver afterwards.

If one forces NDSolve to use a DAE solver by setting SolveDelayed option in the first case, the pdord warning no longer appears, and the resulting "FunctionExpression" won't change:

ndsovletmin = 0.;
ndsolvetmax = 100.;
EoMphi = {phi''[t] + 3. H[t] phi'[t] + phi[t] == 0.};
EoMa = {a'[t] == a[t] H[t], 3. H[t]^2 == phi'[t]^2/2. + phi[t]^2/2., Ne'[t] == H[t]};
bbc = {a[ndsovletmin] == 1., Ne[ndsovletmin] == 0., phi[ndsovletmin] == 2., 
   phi'[ndsovletmin] == 0.};
NDSolve`ProcessEquations[{EoMphi, EoMa, bbc}, {phi[t], phi'[t], a[t], Ne[t]}, {t, 
       ndsovletmin, ndsolvetmax}, SolveDelayed -> #][[1]]["NumericalFunction"][
   "FunctionExpression"] & /@ {Automatic, True}
Function[func, 
  Function @@ {func[[2]] /. Thread[func[[1]] -> (Slot /@ Range@Length@func[[1]])]}] /@ %

(*
{{0. + #5 + 3. #3 #6 + #11, -#2 #3 + #7, 
  3. #3^2 - 0.5 #5^2 - 0.5 #6^2, -#6 + #10, -#3 + #9} &, 
 {0. + #5 + 3. #3 #6 + #11, -#2 #3 + #7, 
  3. #3^2 - 0.5 #5^2 - 0.5 #6^2, -#6 + #10, -#3 + #9} &}
 *)

If one forces NDSolve to use ODE solver, it simply fails:

ndsovletmin = 0.;
ndsolvetmax = 100.;
EoMphi = {ddphi[t] + 3. H[t] dphi[t] + phi[t] == 0., dphi[t] == phi'[t], 
   ddphi[t] == dphi'[t]};
EoMa = {a'[t] == a[t] H[t], 3. H[t]^2 == dphi[t]^2/2. + phi[t]^2/2., Ne'[t] == H[t]};
bbc = {a[ndsovletmin] == 1., Ne[ndsovletmin] == 0., phi[ndsovletmin] == 2., 
   dphi[ndsovletmin] == 0.};
NDSolve[{EoMphi, EoMa, bbc}, {phi[t], phi'[t], a[t], Ne[t]}, {t,
    ndsovletmin, ndsolvetmax}, SolveDelayed -> False]

(* Input returned *)

Also, it's easy to cheat NDSolve in 2nd case. Just add the same derivative term to both sides of the algebraic equations to create 2 fake ODEs:

ndsovletmin = 0.;
ndsolvetmax = 100.;
EoMphi = {dphi[t] == phi'[t], ddphi[t] == dphi'[t]};
EoMa = {a'[t] == a[t] H[t], Ne'[t] == H[t]};
bbc = {a[ndsovletmin] == 1., Ne[ndsovletmin] == 0., phi[ndsovletmin] == 2., 
   dphi[ndsovletmin] == 0.};
fakeODE = {a'[t] + ddphi[t] + 3. H[t] dphi[t] + phi[t] == a'[t] + 0., 
           a'[t] + 3. H[t]^2 == a'[t] + dphi[t]^2/2. + phi[t]^2/2.}
NDSolve[{EoMphi, EoMa, fakeODE, bbc}, {phi[t], phi'[t], a[t], Ne[t]}, {t,
    ndsovletmin, ndsolvetmax}]

NDSolve::pdord

score 1 · Answer 2 · answered Dec 18 '20 at 04:32

Too long for a comment.

I am worried about the ;, I see in your second screenshot. Please do not follow ; with ,

I restart Mathematica so it forgets all previous assignments.

I type literally

ndsovletmin=0.;
ndsolvetmax=100.;
H[t_]:=Sqrt[(phi'[t]^2+phi[t]^2)/6.];
EoMphi={phi''[t]+3. H[t] phi'[t]+phi[t]==0.};
EoMa={a'[t]==a[t] H[t],Ne'[t]==H[t]};
bbc={a[ndsovletmin]==1.,Ne[ndsovletmin]==0.,phi[ndsovletmin]==2.,phi'[ndsovletmin]==0.};
sol=NDSolve[{EoMphi,EoMa,bbc},{phi[t],phi'[t],a[t],Ne[t]},{t,ndsovletmin,ndsolvetmax}]

and it returns

{{phi[t]->InterpolatingFunction[...],phi'[t]->InterpolatingFunction[...],
  a[t]->InterpolatingFunction[...],Ne[t]->InterpolatingFunction[...]}}

with no warnings or errors.

I worry that you solve for phi' instead of only solve for phi and then find the derivative of phi to give you phi' but perhaps your way will not cause serious errors.

What makes a system be differential-algebraic equations?

2 Answers2