# Introduction to Vectors and NDSolve with System of Equations

This question is similar to an older question but with a system of equations.

I have two equations that I want to solve with NDSolve, and here's a similar minimal example:

$$\frac{d\mathbf{F}_1}{d t}=\left[\mathbf{F}_1 \times \mathbf{G}_1 \right]$$ $$\frac{d\mathbf{F}_2}{d t}=\left[\mathbf{F}_2 \times \mathbf{G}_2 \right]$$

$$\mathbf{G}_1 = \mathbf{F}_2 -\hat{i}(\hat{i}\cdot\mathbf{F}_1)$$ $$\mathbf{G}_2 = \mathbf{F}_1 -\hat{i}(\hat{i}\cdot\mathbf{F}_2)$$

Here, $\mathbf{F}_1$ and $\mathbf{F}_1$ are vectors in Cartesian coordinates, $\mathbf{F}_1=F_{1x}\hat{i}+F_{1y}\hat{j}+F_{1z}\hat{k}$ and $\mathbf{F}_2=F_{2x}\hat{i}+F_{2y}\hat{j}+F_{2z}\hat{k}$.

If we code this directly into Mathematica, if gives an error:

g1[t] = f2[t] - {1, 0, 0}*({1, 0, 0}.f1[t])
g2[t] = f1[t] - {1, 0, 0}*({1, 0, 0}.f2[t])
equation = {f1'[t] == Cross[f1[t],g1[t]] , f2'[t] == Cross[f2[t], g2[t]] , f1[0] == {0, 0, -1}, f2[0] == {0, 0, 1}};
sol = NDSolve[equation, {f1, f2}, {t, 0, 10}, MaxSteps -> \[Infinity]]

Cross:The arguments are expected to be vectors of equal length, and the number of arguments is expected to be 1 less than their length.
Cross:The arguments are expected to be vectors of equal length, and the number of arguments is expected to be 1 less than their length.
NDSolve::ndnum: Encountered non-numerical value for a derivative at t == 0..


Perhaps Mathematica is having problems figuring out the dimension of $\mathbf{F}_1$ and $\mathbf{F}_2$? If we remove the last terms with the scalar products,

$$\mathbf{G}_1 = \mathbf{F}_2$$ $$\mathbf{G}_2 = \mathbf{F}_1$$

It solves the code without problem:

g1[t] = f2[t]
g2[t] = f1[t]
equation = {f1'[t] == Cross[f1[t],g1[t]] , f2'[t] == Cross[f2[t], g2[t]] , f1[0] == {0, 0, -1}, f2[0] == {0, 0, 1}};
sol = NDSolve[equation, {f1, f2}, {t, 0, 10}, MaxSteps -> \[Infinity]]

{{f1 -> InterpolatingFunction[{{0., 10.}},


So it works fine without the scalar product term in $\mathbf{G}_1$ and $\mathbf{G}_2$.

What gives? What is causing this problem, and what is correct notation? I would rather do it with vectors than to change it into a huge system of 6 equations.

Your definition for $\mathbf{G}_1$ and $\mathbf{G}_2$ in the code isn't proper. Just observe the output:

g1[t] = f2[t] - {1, 0, 0}*({1, 0, 0}.f1[t]) // Hold // FullForm
(* Hold[Plus[f2[t], Times[-1, Times[List[1, 0, 0], Dot[List[1, 0, 0],
f1[t]]]]]] *)
% // ReleaseHold
(* {-{1, 0, 0}.f1[t] + f2[t], f2[t], f2[t]} *)


As we can see, now f2[t] is everywhere, because Plus has the attribute Listable.

To resolve your problem, we need a function that tries to take the first component of a vector only if it's an explicit vector. This can be achieved by:

Clear@first
first[a_?VectorQ] := {1, 0, 0} a

g1[t] = f2[t] - first@f1[t]
g2[t] = f1[t] - first@f2[t]

equation = {f1'[t] == Cross[f1[t], g1[t]], f2'[t] == Cross[f2[t], g2[t]],
f1[0] == {0, 0, -1}, f2[0] == {0, 0, 1}};

sol = NDSolveValue[equation, {f1, f2}, {t, 0, 10}]


BTW, the following is a more advanced but probably faster way of defining first:

first = Compile[{{a, _Real, 1}}, {1, 0, 0} a,
RuntimeOptions -> EvaluateSymbolically -> False]

• This is great. It seems like first[a_?VectorQ] := {1, 0, 0} a can be changed to {0, 1, 0} or whatever direction as well. Thanks for this. Jul 3, 2018 at 14:01
• How do I extract the solution from sol? f1[t] /. sol[[1]]` doesn't seem to work. Should I ask this as another question? Jul 7, 2018 at 13:15
• @axsvl77 I can edit my answer to include a solution, but my solution is not perfect, and AFAIK this non-trivial question hasn't been asked in this site before, so I suggest starting a new question. Jul 7, 2018 at 13:23
• I have added a question to address this issue Jul 7, 2018 at 13:50