# Vector-valued ParametricNDSolve, solving for a combination

I am integrating two vector differential equations using ParametricNDSolve, one for $$\mathbf{Y}$$ and one for $$\mathbf{Z}$$, and then I'm interested in a combination of the two of them. ParametricNDSolve allows for specifying the output to just be a function of the two, for instance $$\mathbf{Y} \cdot \mathbf{Z}$$, but the combination I want is (for the case of sixth order vector equations): $$\psi = y_1 z_6 -y_2 z_5 + y_3 z_4 + y_4 z_3 - y_5 z_2 +y_6 z_1$$

Is there a way to get ParametricNDSolve to give me this combination directly?

Example code:

A = {{0, 1, 0, 0}, {0, 0, 1, 0}, {0, 0, 0, 1}, {q, 0, -Sin[x], 0}};
Q[x_, q_] = {{0, 1, 0, 0, 0, 0}, {0, 0, 1, 1, 0, 0}, {0, -Sin[x],0, 0, 1, 0},
{0, 0, 0, 0, 1, 0}, {-q, 0, 0, -Sin[x], 0, 1}, {0, -q, 0, 0, 0, 0}};
n = Length[Q[x, q]];
cA = With[{code = N@A}, Compile[{{x, _Real}, {q, _Real}}, code]];
F[x_?NumericQ, q_?NumericQ, Y_?VectorQ] :=
(Q[x, q] - DiagonalMatrix[ConstantArray[1, n]]).Y;
xa = -4; xb = 4; xm = 0;
k = -((-xb + xm)/(xa - xm)); c = -((-xa xm + xb xm)/(xa - xm));
Y0 = {0, 0, 0, 0, -1, 0};
Z0 = {0, 1, 0, 0, 0, 0};

YZsol = ParametricNDSolveValue[{Y'[x] == F[x, q, Y[x]], Y[xa] == Y0,
Z'[x] == 1/k F[k x + c , q, Z[x]], Z[xa] == Z0}, {Y, Z}, {x, xa, xm}, q];


I'm interested in finding where $$\psi=0$$, which I can do with this construction, but is there a more direct way to get it out of the ParametricNDSolveValue?

psi[q_?NumericQ] := (Y Z - Y Z + Y Z + Y Z -
Y Z + Y Z) /. {Y[a_] :> YZsol[q][][xm][[a]],
Z[a_] :> YZsol[q][][xm][[a]]};

FindRoot[psi[q], {q, 3}]


Edit:

Carl's answer does exactly what I need, but I see a strange timing effect, where after I use FindRoot it takes 6 times longer to do the same calculation:

Clear[YZsol2]; YZsol2 = ParametricNDSolveValue[{Y'[x] == F[x, q, Y[x]], Y[xa] == Y0, Z'[x] == 1/k F[k x + c, q, Z[x]], Z[xa] == Z0}, Y[xm].form.Z[xm], {x, xa, xm}, q];
AbsoluteTiming[YZsol2 /@ Range]
(* {0.259752, Null} *)

Clear[YZsol3]; YZsol3 =
ParametricNDSolveValue[{Y'[x] == F[x, q, Y[x]], Y[xa] == Y0,
Z'[x] == 1/k F[k x + c, q, Z[x]], Z[xa] == Z0}, Y[xm].form.Z[xm], {x, xa, xm}, q];
FindRoot[YZsol3[q], {q, 3}];
AbsoluteTiming[YZsol3 /@ Range]
(* {1.69719, Null} *)


## 2 Answers

Define a quadratic form matrix:

form = {
{0, 0, 0, 0, 0, 1},
{0, 0, 0, 0, -1, 0},
{0, 0, 0, 1, 0, 0},
{0, 0, 1, 0, 0, 0},
{0, -1, 0, 0, 0, 0},
{1, 0, 0, 0, 0, 0}
};


Then:

YZsol2 = ParametricNDSolveValue[
{Y'[x] == F[x,q,Y[x]], Y[xa]==Y0,
Z'[x] == 1/k F[k x+c,q,Z[x]], Z[xa]==Z0},
Y[xm] . form . Z[xm],
{x,xa,xm},
q
];


will produce the desired function. For example:

YZsol2 /@ Range[0, 4]
psi /@ Range[0, 4]


{-24.1853, -148.784, -6240.84, -7841.06, 12815.}

{-24.1853, -148.784, -6240.84, -7841.06, 12815.}

Another possibility is to include your desired output as an additional equation (constraint):

YZsol3 = ParametricNDSolveValue[
{
Y'[x] == F[x,q,Y[x]], Y[xa]==Y0,
Z'[x] == 1/k F[k x+c,q,Z[x]], Z[xa]==Z0,
U[x] == Y[x] . form . Z[x]
},
U[xm],
{x,xa,xm},
q
];


Check:

YZsol3 /@ Range[0, 4]


{-24.1853, -148.784, -6240.84, -7841.06, 12815.}

Note that there is a bug as of June 2019 with using this inside a FindRoot, which makes it much slower unless we add Method -> {"ParametricSensitivity" -> None} as detailed in Vector ParametricNDSolve and FindRoot interaction).

• Thanks, of course that is a way to write it as a vector/matrix equation – KraZug May 31 at 14:31
• Huh, that works really well, except they are both much slower than the original formulation. The constraint version also fails to find the root (presumably because it is using a different internal method on the NDSolve. – KraZug May 31 at 14:42
• I edited your answer to include a link to another question where I show the workaround to not be slow with FindRoot. – KraZug Jun 7 at 7:35

A hybrid of yours and @CarlWoll's approach seems to perform well after FindRoot.

Clear[YZsol, psi2]
form = Reverse@DiagonalMatrix@PadRight[{}, 6, {1, -1, 1}];
YZsol = ParametricNDSolveValue[{Y'[x] == F[x, q, Y[x]], Y[xa] == Y0,
Z'[x] == 1/k F[k x + c, q, Z[x]], Z[xa] == Z0}, {Y, Z}, {x, xa,
xm}, q];
psi2[q_?NumericQ] := (Y.form.Z) /. {Y :> YZsol[q][][xm],
Z :> YZsol[q][][xm]};
FindRoot[psi2[q], {q, 3}]
AbsoluteTiming[psi2 /@ Range]

Clear[YZsol2];
YZsol2 = ParametricNDSolveValue[{Y'[x] == F[x, q, Y[x]], Y[xa] == Y0,
Z'[x] == 1/k F[k x + c, q, Z[x]], Z[xa] == Z0},
Y[xm].form.Z[xm], {x, xa, xm}, q];
AbsoluteTiming[YZsol2 /@ Range]

Clear[YZsol3]; YZsol3 =
ParametricNDSolveValue[{Y'[x] == F[x, q, Y[x]], Y[xa] == Y0,
Z'[x] == 1/k F[k x + c, q, Z[x]], Z[xa] == Z0},
Y[xm].form.Z[xm], {x, xa, xm}, q];
FindRoot[YZsol3[q], {q, 3}];
AbsoluteTiming[YZsol3 /@ Range]
(* {q -> 3.55393} *)
(* {0.279002, {-148.784, -6240.84, -7841.06, 12815., 63393.7,
134180., 197976., 214749.,
139054., -71335.6, -446840., -998422., -1.71155*10^6, -2.5422*10^6, \
-3.41516*10^6, -4.22468*10^6, -4.83742*10^6, -5.09767*10^6, \
-4.83432*10^6, -3.86961*10^6}} *)
(* {0.337883, {-148.784, -6240.84, -7841.06, 12815., 63393.7,
134180., 197976., 214749.,
139054., -71335.6, -446840., -998422., -1.71155*10^6, -2.5422*10^6, \
-3.41516*10^6, -4.22468*10^6, -4.83742*10^6, -5.09767*10^6, \
-4.83432*10^6, -3.86961*10^6}} *)
(* {2.20664, {-148.785, -6240.83, -7841.05, 12815., 63393.6,
134180., 197976., 214748.,
139054., -71335.6, -446840., -998421., -1.71155*10^6, -2.5422*10^6, \
-3.41516*10^6, -4.22467*10^6, -4.83742*10^6, -5.09767*10^6, \
-4.83432*10^6, -3.86961*10^6}} *)

• Er… what's the difference between YZsol2 and YZsol3? – xzczd Jun 1 at 4:53
• @xzczd, there is no difference, just one has FindRoot applied first – KraZug Jun 1 at 6:40
• @KraZug So, if FindRoot is called first, the corresponding ParametricFunction will be slower? – xzczd Jun 1 at 9:09
• This is basically the same as my approach directly. I was thinking that using ParametricNDSolve to calculate the quantity of interest would be faster and perhaps more accurate. – KraZug Jun 1 at 9:11
• @xzczd, exactly. And the FindRoot is slower too. It is marginally faster to evaluate individual points until you use FindRoot , then it slows down considerably – KraZug Jun 1 at 9:13