2
$\begingroup$

I'm trying to solve this system of ODEs numerically:

$$\vec{Q} = \vec{P} \times \vec{Q}$$ where vectors $\vec{P}, \vec{Q}$ have 3 components:

$$\vec{Q} = \left( Q_1, Q_2, Q_3 \right) $$ $$\vec{P} = \left( P_1, P_2, P_3 \right)$$ with given $$P_1 = 1$$ $$P_2 = a \cos \left[w(-\log(|t|))\right] + b \sin \left[w(-\log(|t|))\right]$$ $$P_3 = -a \sin \left[w(-\log(|t|))\right] + b \cos \left[w(-\log(|t|))\right]$$ Because of rotation character of this equation $$|\vec{Q}| = const$$

But when I use NDSolve for this and plot the result, I obtain: enter image description here

As we see, $|\vec{Q}|$ is not constant (it's growing with time).

So, I think that the reason is accumulating errors during calculations. Here's a code:

 sol = NDSolve[{Q1'[t] == ((a*Cos[w*-Log[Abs[t]]] + b*Sin[w*-Log[Abs[t]]])*
   Q3[t] + (a*Sin[w*-Log[Abs[t]]] - b*Cos[w*-Log[Abs[t]]])*Q2[t]), 
 Q2'[t] == ((b*Cos[w*-Log[Abs[t]]] + a*Sin[w*-Log[Abs[t]]])*Q1[t] - Q3[t]), 
 Q3'[t] == ((b*Sin[w*-Log[Abs[t]]] + a*Cos[w*-Log[Abs[t]]])*Q1[t] + Q2[t]),
 Q1[t1] == 0, Q2[t1] == 0, Q3[t1] == 1}, {Q1[t], Q2[t], Q3[t]}, {t, t1, t3},
 AccuracyGoal -> 20, PrecisionGoal -> 20, WorkingPrecision -> 35, MaxSteps->
 Infinity] 

 Plot[{Evaluate[Q1[t]] /. sol, Evaluate[Q2[t]] /. sol, Evaluate[Q3[t]] /.
 sol},  {t, t1, t3}, PlotRange -> 1.5, AxesLabel -> {t, Q3[t]}]

I even wrote all used values ($a$, $b$, etc.) with explicit precision (like a = 0.5`20).

If I increase AccuracyGoal and PrecisionGoal (both to value 30 and WorkingPrecision 35) I get something very wrong:

enter image description here

and I don't know why this occurs.

UPD: typical value of $w$ is about $20.0$.

UPD2: I copied Winther's code and got errors:

 NDSolve::underdet: 
 There are more dependent variables, {Q1[t],Q2[t],Q3[t]}, than \
 equations, so the system is underdetermined. 
 NDSolve::dsvar: 0.006021326428571428` cannot be used as a variable. 
 ReplaceAll::reps: "{NDSolve[<<1>>]} is neither a list of replacement rules
 nor a valid dispatch table, and so cannot be used for replacing."
 NDSolve::dsvar: 0.006021326428571428` cannot be used as a variable.
 ReplaceAll::reps: "{NDSolve[<<1>>]} is neither a list of replacement rule
 nor a valid dispatch table, and so cannot be used for replacing."
 NDSolve::dsvar: 1.026327448877551` cannot be used as a variable. 

and so on.

All values have been cleared and initialized again before compilation.

$\endgroup$
  • $\begingroup$ What's w supposed to be? $\endgroup$ – J. M. is away Jun 29 '15 at 7:43
  • $\begingroup$ w is a ''frequency'' of rotation of vector P $\endgroup$ – newt Jun 29 '15 at 8:03
  • 2
    $\begingroup$ Okay, I was hinting that you include a typical value of w in the question… $\endgroup$ – J. M. is away Jun 29 '15 at 8:06
  • $\begingroup$ added, typical value is about 20.0 $\endgroup$ – newt Jun 29 '15 at 10:41
  • $\begingroup$ I would first check how well the solution satisfies the ODE. I would also tend to keep AccuracyGoal and PrecisionGoal constant while increasing WorkingPrecision; then one might expect round-off error to diminish. $\endgroup$ – Michael E2 Jun 29 '15 at 18:45
4
$\begingroup$

I rewrote your code and it works for me so I suspect there is a bug in your code.

Some tips: when you write out the equation as you have done it is very easy to make mistakes. What I recomment is that you define $P=\{P_1,P_2,P_3\}$ and $Q=\{Q_1,Q_2,Q_3\}$ only once (instead of retyping it every time you need it) and then write the equation directly in terms of these vectors using Mathematica's vector methods. This can be done quite easily as shown below:

 (* Define some constants *)
 t1 = 0.5; t3 = 20; w = 20.0; a = 0.5; b = 20.0;

 (* Define P *)
 P1 = 1;
 P2 =  a Cos[z] + b Sin[z];
 P3 = -a Sin[z] + b Cos[z];
 P = {P1, P2, P3} /. z -> -w Log[Abs[t]];

 (* Define Q *)
 Q = {Q1[t], Q2[t], Q3[t]};

 (* Define equation(s) *)
 eq = D[Q,t] - Cross[P,Q];
 Do[eq[[i]] = eq[[i]] == 0, {i, 1, Length[eq]}];

 (* Define initial condition(s) *)
 ic = {0, 0, 1};
 Do[ic[[i]] = (Q /. t -> t1)[[i]] == ic[[i]], {i, 1, Length[ic]}]

 (* Solve equation(s) *)
 sol = NDSolve[{eq, ic}, Q, {t, t1, t3}];

 Plot[{Norm[Q] /. sol}, {t, t1, t3}, 
    PlotRange -> 1.5, AxesLabel -> {t, "|Q[t]|"}]

Below is a plot of $|Q(t)|$ for a run of the code above. As you can see it works just as desired:

$~~~~~~~~~~~~~~~~~~~~~~~~~~~$enter image description here

$\endgroup$
  • $\begingroup$ Thanks. But I copied your code and got error (see updated question) $\endgroup$ – newt Jul 3 '15 at 12:32
  • $\begingroup$ @newt I made a mistake in writing in the code. The code I have now works when I copy it into my Mathematica 9 $\endgroup$ – Winther Jul 3 '15 at 12:52
  • $\begingroup$ Ah, thanks. That's how it works. $\endgroup$ – newt Jul 3 '15 at 13:23

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.