0
$\begingroup$

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.

Can anyone please help to explain why this is occurring? Thanks a lot.

$\endgroup$
2
  • 2
    $\begingroup$ 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. $\endgroup$
    – Nasser
    Commented Feb 6, 2023 at 9:14
  • $\begingroup$ @Nasser Thank you for the clarification. So it is safe to say NDSolve can't handle Differential-algebraic BVPs directly. $\endgroup$
    – Jack L
    Commented Feb 6, 2023 at 12:26

1 Answer 1

0
$\begingroup$

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.} *)

enter image description here

$\endgroup$
1
  • $\begingroup$ Thanks a lot for your guidance! It works like a charm. $\endgroup$
    – Jack L
    Commented Feb 6, 2023 at 12:23

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.