4
$\begingroup$

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[1] Z[6] - Y[2] Z[5] + Y[3] Z[4] + Y[4] Z[3] - 
     Y[5] Z[2] + Y[6] Z[1]) /. {Y[a_] :> YZsol[q][[1]][xm][[a]], 
    Z[a_] :> YZsol[q][[2]][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[20]] 
(* {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[20]]
(* {1.69719, Null} *)
$\endgroup$
5
$\begingroup$

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

$\endgroup$
  • $\begingroup$ Thanks, of course that is a way to write it as a vector/matrix equation $\endgroup$ – KraZug May 31 at 14:31
  • $\begingroup$ 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. $\endgroup$ – KraZug May 31 at 14:42
  • $\begingroup$ I edited your answer to include a link to another question where I show the workaround to not be slow with FindRoot. $\endgroup$ – KraZug Jun 7 at 7:35
2
$\begingroup$

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][[1]][xm], 
    Z :> YZsol[q][[2]][xm]};
FindRoot[psi2[q], {q, 3}] 
AbsoluteTiming[psi2 /@ Range[20]]

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[20]]


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[20]]
(* {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}} *)
$\endgroup$
  • $\begingroup$ Er… what's the difference between YZsol2 and YZsol3? $\endgroup$ – xzczd Jun 1 at 4:53
  • $\begingroup$ @xzczd, there is no difference, just one has FindRoot applied first $\endgroup$ – KraZug Jun 1 at 6:40
  • $\begingroup$ @KraZug So, if FindRoot is called first, the corresponding ParametricFunction will be slower? $\endgroup$ – xzczd Jun 1 at 9:09
  • $\begingroup$ 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. $\endgroup$ – KraZug Jun 1 at 9:11
  • $\begingroup$ @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 $\endgroup$ – KraZug Jun 1 at 9:13

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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