4
$\begingroup$

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.

$\endgroup$
4
$\begingroup$

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]
$\endgroup$
  • $\begingroup$ 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. $\endgroup$ – axsvl77 Jul 3 '18 at 14:01
  • $\begingroup$ How do I extract the solution from sol? f1[t] /. sol[[1]] doesn't seem to work. Should I ask this as another question? $\endgroup$ – axsvl77 Jul 7 '18 at 13:15
  • 1
    $\begingroup$ @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. $\endgroup$ – xzczd Jul 7 '18 at 13:23
  • $\begingroup$ I have added a question to address this issue $\endgroup$ – axsvl77 Jul 7 '18 at 13:50

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.