5

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.

xzczd
  • 65,995
  • 9
  • 163
  • 468
Astor Florida
  • 1,326
  • 7
  • 21

1 Answers1

5

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]
xzczd
  • 65,995
  • 9
  • 163
  • 468