NDSolve mixing many scalar and vector equations

Question

I have a set of scalar equations in many unknowns, which I want to combine with a vector equation inside an NDSolve. The equations are a mix of differential and algebraic equations.

A set of scalar equations:

eqns = {-a[0][t] + a[1][t] == 0, a[1]'[t] == - fx[1][t] + fy[1][t] - a[1][t],
        a[2]'[t] == - fx[2][t] + fy[2][t] - a[2][t], a[2][t] + a[3][t] == 0,
        0 == x[0][t], a[1][t] == - x[0][t] + x[2][t], a[2][t] == - x[1][t] + x[3][t], 
        a[3][t] == x[1][t] + x[3][t], 1 == y[0][t], a[1][t] == - y[0][t] + y[2][t], 
        a[2][t] == - y[1][t] + y[3][t], a[3][t] == -y[2][t] + y[3][t]};

Along with a vector equation, defined using a function fCalc.

feqn = {{fx[0][t], fy[0][t]}, {fx[1][t], fy[1][t]}, {fx[2][t], fy[2][t]}, {fx[3][t], fy[3][t]}} == 
      fCalc[{{x[0][t], y[0][t]}, {x[1][t], y[1][t]}, {x[2][t], y[2][t]}, {x[3][t], y[3][t]}}];

My actual function fCalc is complicated, but it only evaluates for numerical input. It takes a list of {x,y} points and returns a list of {x,y} points, for instance:

fCalc[pts_ /; MatrixQ[pts, NumericQ]] := pts

Some initial conditions and the variables used:

initcs = {a[0][0] == 0, a[1][0] == 0, a[2][0] == 0, a[3][0] == 0};
vars=Flatten[Table[{a[j], x[j], y[j], fx[j], fy[j]}, {j, 0, 3}]];

Trying to solve this directly gives an error, as Mathematica tries to evaluate this vector function before starting, and refuses to start as it doesn't see enough equations for the unknowns.

NDSolve[{eqns, feqn, initcs},vars, {t, 0, 1}];
NDSolve::underdet: There are more dependent variables, than equations, so the system is underdetermined.

Removing the numerical requirement the system is well-defined:

fCalc2[pts_] := pts
AbsoluteTiming[NDSolve[{eqns, feqn /. fCalc -> fCalc2, initcs}, vars, {t, 0, 1}];]

So my question is, how to I trick Mathematica into waiting to evaluate the NDSolve until after giving values to the variables, and without having to repeatedly evaluate the fCalc function. There is an answer for solving for a vector system at this question, but I can't write my whole system in terms of a derivative of a vector as in that case.

My actual system is much more complicated (see the initial version of this question if you want to see), but I definitely need the numerical restriction on fCalc.

Answer 1

The system of equations can be solved following the procedure in the accepted answer to question 78641. First, observe that the ODEs in eqns all have the form,

a[i]'[t] == - fx[i][t] + fy[i][t] - a[i][t]

and so can be written in vector form,

a'[t] == - fx[t] + fy[t] - a[t]

Consequently, the solution is given by

NDSolve[{a'[t] == -a[t] - fCalc[a[t]], a[0] == {0, 0}}, a, {t, 0, 1}]

where the function fCalc[a[t]] returns the array of values fx[t] - fy[t] corresponding to a[t]. This function has the form,

fCalc[pts : {_?NumberQ ..}] := Module[{l = Length[pts], s}, Do[a[i] = pts[[i]], {i, l}]; 
    s = Solve[eqsim, varsim] // Flatten; ...]

where ... represents whatever code the OP uses to compute fx - fy, given values of x and y. For instance, if the process were simply to set fx equal to x and fy equal to y, then ... would be Table[(x[i] - y[i]) /. s, {i, l}]. eqsim is eqns with the ODEs removed, and varsim is a List of the x and y variables and the a variables not appearing in the ODEs.

For the equations given in the question, this reduces to

Clear[a]
eqsim = Delete[eqns, {{2}, {3}}] /. {a[i_][t] -> a[i], x[i_][t] -> x[i], y[i_][t] -> y[i]};
varsim = Join[Flatten[Table[{x[j], y[j]}, {j, 0, 3}]], {a[0], a[3]}];
fCalc[pts : {_?NumberQ ..}] := Module[{l = Length[pts], s}, Do[a[i] = pts[[i]], {i, l}]; 
    s = Solve[eqsim, varsim] // Flatten; Table[(x[i] - y[i]) /. s, {i, l}]];

NDSolve[{a'[t] == -a[t] - fCalc[a[t]], a[0] == {0, 0}}, a, {t, 0, 1}];
Plot[a[t] /. sol, {t, 0, 1}, AxesLabel -> {t, a}]

(The curves for a[1] and a[2] coincide.)