# Error when using NDsolve for Differential-Algebraic Equations

For the system of differential-algebraic equations with boundary conditions,

eqns = {D[p[x], x, x] - q[x] == Sin[x], p[x] + q[x] == 1};
bcs = {p[1/2] == 0, D[p[x], x] == 0 /. x -> 0};
sol3 = NDSolve[{eqns, bcs}, {p[x], q[x]}, x]


Mathematica returns an error:

NDSolve::bvdae: Differential-algebraic equations must be given as initial value problems.


• the error seems to be clear. It says Differential-algebraic equations must be given as initial value problems. but you have given {p[1/2] == 0, D[p[x], x] == 0 /. x -> 0} i.e. both initial conditions are not at the same point. This is considered not an initial value problem but a BVP. an IVP problem will have all the conditions given at same location. Commented Feb 6, 2023 at 9:14
• @Nasser Thank you for the clarification. So it is safe to say NDSolve can't handle Differential-algebraic BVPs directly. Commented Feb 6, 2023 at 12:26

You may us the "shooting method" to find an initial value for p, so that p[1/2]==0. Toward this aim, we define a function "fun" that takes an initial value and returns the value p[1/2] corresponding to the given inital value. Then we can use "FindRoot" to find the correct inital value.

Clear["Global*"]

fun[y_?NumericQ] := Module[{eqns, sol, ics, x, p, q},
eqns = {p''[x] - q[x] == Sin[x], p[x] + q[x] == 1};
ics = {p[0] == y, p'[x] == 0 /. x -> 0};
sol = NDSolve[{eqns, ics}, {p[x], q[x]}, {x, 0, 1/2}][[1]];
p[x] /. sol /. x -> 1/2];

FindRoot[fun[x], {x, 0}]

(*{x -> -0.162645}*)


To test we can calculate the solutions for t=0..10 and check if p[1/2] is zero. (Note NDSolve needs an interval for the independent variable):

eqns = {p''[x] - q[x] == Sin[x], p[x] + q[x] == 1};
ics = {p[0] == -0.162645113335746, p'[0] == 0};
sol[x_] = {p[x], q[x]} /.
NDSolve[{eqns, ics}, {p[x], q[x]}, {x, 0, 10}][[1]] ;

sol[1/2]
Plot[Evaluate[sol[x]], {x, 0, 10}]

(* {4.91608*10^-8, 1.} *)


• Thanks a lot for your guidance! It works like a charm. Commented Feb 6, 2023 at 12:23