3

I am attempting to solve a system of ODEs where a, b, cv1 and cv2 are vectors.

Tu = {{-3 / 4, 1 / 4}, {5, -6}}`
tu = {2/4, 1}
pu = {.7, .3}

a[y_] = {Subscript[a, 1][y], Subscript[a, 2][y]};
b[y_] = {Subscript[b, 1][y], Subscript[b, 2][y]};
cv1[y_] = {Subscript[cv1, 1][y], Subscript[cv1, 2][y]}
cv2[y_] = {Subscript[Subscript[cv2, 1][y], 1], Subscript[Subscript[cv2, 2][y], 2]}


(* Ode system )
a'[y] == a[y].Tu
b'[y] == Tu.b[y]
cv1'[y] == Tu.cv1[y] + Subscript[a, 1][y].tu
cv2'[y] == Tu.cv2[y] + Subscript[a, 2][y].tu)


S = 
  NDSolve[
    {a'[y] == a[y].Tu,
     b'[y] == Tu.b[y],
     cv1'[y] == Tu.cv1[y] + Subscript[a, 1][y].tu,
     cv2'[y] == Tu.cv2[y] + Subscript[a, 2][y].tu, a[0] == pu, 
     b[0] == tu, cv1[0] == cv2[0] == 10^-5}, 
   {a[y], b[y], cv1[y], cv2[y] }, {y, 0, ∞}]

I am getting an output informing me that there (incorrectly) exist more dependent variables than equations. How should I specify this system?

Edit

I fixed the initial code (the above is now the updated version); the previous error no longer occurs, though I now receive the following:

NDSolve::ndnco: The number of constraints (6) (initial conditions) is not equal to the total differential order of the system plus the number of discrete variables (8).

EDIT:

Note that using the system as specified in the answer will not retrieve the correct cv1 and cv2 expressions. Instead specify those in the system as:

cv1'[y] == 
 Tu.cv1[y] + (MapIndexed[Indexed[a[y], #2] &, Range[2]] // First)*tu,
cv2'[y] == 
 Tu.cv2[y] + (MapIndexed[Indexed[a[y], #2] &, Range[2]] // Last)*tu
Community
  • 1
BayesIsBaye
  • 145
  • 5
  • There is no initial data. – Alex Trounev Jan 22 '19 at 00:01
  • My bad, I updated it with the correct initial data and am getting new errors as will be listed – BayesIsBaye Jan 22 '19 at 01:04
  • It's just unclear to me what system you're trying to solve. Is Subscript[a, 1][y] and Subscript[a, 2][y] scalar? If so, why there exists Subscript[a, 1][y].tu in the system? (Dot (.) is for products of vectors, matrices, and tensors. ) Is cv2 just a vector like a, b and cv1? If so, why is it defined in such a strange way? Can you show us the system in traditional math notation? – xzczd Jan 25 '19 at 05:26

1 Answers1

3
Tu = {{-3/4, 1/4}, {5, -6}};

tu = {2/4, 1};
pu = {.7, .3};
ym = 1;


(*Ode system*)
eq = {a'[y] == a[y].Tu, b'[y] == Tu.b[y], 
   cv1'[y] == Tu.cv1[y] + First[a[y]]*tu, 
   cv2'[y] == Tu.cv2[y] + Last[a[y]]*tu};
ic = {a[0] == pu, b[0] == tu, cv1[0] == cv2[0] == {10^-5, 10^-5}};



S = NDSolve[{eq, ic}, {a, b, cv1, cv2}, {y, 0, ym}]

{Plot[a[y] /. S, {y, 0, ym}, AxesLabel -> {"y", "a"}], 
 Plot[b[y] /. S, {y, 0, ym}, AxesLabel -> {"y", "b"}, 
  PlotRange -> All], 
 Plot[cv1[y] /. S, {y, 0, ym}, AxesLabel -> {"y", "cv1"}], 
 Plot[cv2[y] /. S, {y, 0, ym}, AxesLabel -> {"y", "cv2"}]}

fig1

If we change the code, excluding First and Last, the result will change, although in symbolic form the equations look the same.

(*Ode system*)
eq = {a'[y] == a[y].Tu, b'[y] == Tu.b[y], 
   cv1'[y] == Tu.cv1[y] + (a[y].{1, 0})*tu, 
   cv2'[y] == Tu.cv2[y] + (a[y].{0, 1})*tu};
ic = {a[0] == pu, b[0] == tu, cv1[0] == cv2[0] == {10^-5, 10^-5}};
S = NDSolve[{eq, ic}, {a, b, cv1, cv2}, {y, 0, ym}];
{Plot[a[y] /. S, {y, 0, ym}, AxesLabel -> {"y", "a"}], 
 Plot[b[y] /. S, {y, 0, ym}, AxesLabel -> {"y", "b"}, 
  PlotRange -> All], 
 Plot[cv1[y] /. S, {y, 0, ym}, AxesLabel -> {"y", "cv1"}], 
 Plot[cv2[y] /. S, {y, 0, ym}, AxesLabel -> {"y", "cv2"}]}

fig2

I assume that functions First and Last are ignored when solving equations. @Janeiro offers its solution.

Alex Trounev
  • 44,369
  • 3
  • 48
  • 106
  • A colleague who had coded this in R had noted that he had different forms for CV1 and CV2 -- is the First[] and Last[] commands here proper in the system of equations wrt a's elements? – BayesIsBaye Jan 24 '19 at 16:57
  • 1
    Problem solved. Use : cv1'[y] == Tu.cv1[y] + (MapIndexed[Indexed[a[y], #2] &, Range[2]] // First)tu, cv2'[y] == Tu.cv2[y] + (MapIndexed[Indexed[a[y], #2] &, Range[2]] // Last)tu (will update question) – BayesIsBaye Jan 24 '19 at 17:14
  • The result is the same. What is the point of adding? – Alex Trounev Jan 24 '19 at 18:33
  • @Janeiro You are right, really NDSolve[] ignores First[] and Last[]. I published the second solution. – Alex Trounev Jan 24 '19 at 19:20
  • NDSolve doesn't ignore First and Last, because they've already disappeared when the eq=… line is executed: First@a[y] evaluates to y immediately. Notice most functions for list manipulation are not limited to List manipulation e.g. try {First@f[x], Most@f[x,y,z]}. – xzczd Jan 25 '19 at 05:08
  • @xzczd We can eliminate eq and put the equations right in NDSolve[]. Then we will see that NDSolve ignores that First[] and Last[] acts on the list a[y] and takes a[y] as a scalar. Although in calculations a[y] is taken as a vector. This is ridiculous. – Alex Trounev Jan 25 '19 at 12:01
  • Eliminating eq doesn't change anything, because NDSolve doesn't have attribute HoldAll or HoldFirst etc. so the first argument of NDSolve i.e. the equation system will be evaluated before sending into NDSolve. Related: https://mathematica.stackexchange.com/a/97006/1871 – xzczd Jan 25 '19 at 15:17
  • @xzczd Thanks, I understand, but it is funny that the function in the same row is considered both as a vector and as a scalar. – Alex Trounev Jan 25 '19 at 16:35