# NDSolve mixing many scalar and vector equations

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.

• Can you provide a simplified version of your code that still illustrates the essential issues? Few readers will undertake to address so much code. – bbgodfrey Apr 18 '16 at 22:58
• @bbgodfrey, I have condensed it further, I'm loathe to go further as I think it'll lose the actual structure that I need. – KraZug Apr 19 '16 at 13:46
• You can try a numerical solution like this – Sumit May 12 '16 at 8:56
• Yes, I could do the timestepping myself, but I would like to use NDSolve so I don't have to – KraZug May 15 '16 at 4:56
• Hmm, nobody have an answer for me? There is a sweet sweet bounty available... – KraZug May 16 '16 at 8:56

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

• I thought of this, but I was searching for a way to get NDSolve to handle mixed dimensions. – Michael E2 May 18 '16 at 0:42
• @MichaelE2 If you cannot find a better way, I do not know who could. Certainly, I tried, but without success. – bbgodfrey May 18 '16 at 2:18
• Thanks for having a go. I think I may be able to adapt this to my case, with having all the algebraic equations solved inside the vector function. In general, I have a'[t] == (finite difference approximation of a spatial derivative, including pieces of a,x,y), not having the nice structure that the reduced example did. – KraZug May 18 '16 at 18:43
• @KraZug The entire right side of the ODE for a[t] can be computed by fCalc[], so the method in my answer should work in general. I agree that a more elegant approach would be preferable, but I (and Michael E2, it would appear) each tried hard to find one but without success. – bbgodfrey May 18 '16 at 18:54
• I'll give it a go next week, I'm away for a few days. Thanks for your help (and @MichaelE2). – KraZug May 19 '16 at 6:37