Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

I'm trying to set up a system of differential equations for passing to NDSolve. Note that my initial conditions are vector valued so Mathematica should know that raptor[[i]] is a vector valued function. However, Mathematica seems to think that that raptor[[i]] is a scalar valued function (this is indicated by my output, which is included below).

Input:

human[\[Theta]_, t_] := 0 t*{Cos[\[Theta]], Sin[\[Theta]]} + {0, 10 Sin[\[Pi]/3]}

V = {25, 10, 25};
P = {{-10, 0}, {0, 20 Sin[\[Pi]/3]}, {10, 0}};
\[Theta] = 0;

 Flatten[
  Table[
  {
   raptor[[i]]'[t] == Normalize[human[\[Theta], t] - raptor[[i]][t]]
   raptor[[i]][0] == P[[i]]
  }, {i, 3}]
 ]

Output:

{ Raptor[[1]]'[t]=={ first component of human - raptor[[i]][t],
                     second component of human - raptor[[i]][t]},
  Raptor[[1]][0] =={-10,0}, 
  Raptor[[2]]'[t]== ... 
}

Furthermore, changing the raptor indexes into subscripts does not change anything. When I try to use this system in NDSolve I get the error "Computed derivatives do not have dimensionality consistent with the initial conditions."

After talking to my TA for the course, we got this to work by defining human[theta, t] as one of the differential equations. However, my TA didn't know why my code failed (it seemed very reasonable to both of us) and said I could ask the Mathematica stack exchange.

So, could anyone explain why my code fails? Namely, why does Mathematica think that raptor[[i]] is a scalar valued function even though its initial condition is clearly vector valued? Is there a better way to do this then putting human[theta, t] into the list of equations for NDSolve?

PS: for those wondering, this is solving a variant of problem #2 from http://xkcd.com/135/

share|improve this question
    
+1 for velociraptor –  Silvia Apr 23 '13 at 3:54
    
If you don't consider the acceleration, the bottom raptors will get you before you can take any advantage of the top raptor's broken leg. Sorry. –  belisarius Apr 23 '13 at 3:57
    
I don't think running in straight line is the best strategy. But well, it only ask to "maximize the time you stay alive", after that, ... –  Silvia Apr 23 '13 at 3:59
1  
I presume you've looked up "pursuit curves"? –  J. M. Apr 23 '13 at 3:59
1  
I think you should initialize raptor as something like raptor = {raptor1, raptor2, raptor3} before using raptor[[i]]. –  Silvia Apr 23 '13 at 4:11

2 Answers 2

The following example is written for Mathematica 9, and will NOT be compatible with previous versions.

According to your question, this looks like some kind of homework, so I'd show a slightly different example, just to demonstrate how to index variables (by one of many possible ways), and how to use some features of NDSolve that may be helpful for chasing problems.

So here we have a target initially at $(0, 0)$, and some chasers distribute uniformly on a circle centered at $(0,0)$ with radius $10$:

Clear[targetPos]
targetPos[v_, θ_, t_] := v t*{Cos[θ], Sin[θ]}

chaserNumber = 6;

chaserOrigPos = 
  10 {Cos[#], Sin[#]} & /@ Range[0, 2 π - (2 π)/chaserNumber, (2 π)/chaserNumber];

Then we define the ODE system. Please notice how chaserPos is indexed. For something like $f_{i,j}(x)$, instead of using f[[i,j]][x] or even Subscript[f,i,j][x], f[i][j][x] might be a more convenient choice (you don't need to initialize it and don't have upper/lower limits for indices):

Clear[eqsys]
eqsys[chaserVelocity_, targetVelocity_, θ_] := Flatten[Table[{
    Derivative[1][chaserPos[i][j]][
      t] == ((chaserVelocity #)/Sqrt[#.#] &)[
       targetPos[targetVelocity, θ, t] - Through[(chaserPos[i] /@ {1, 2})[t]]
       ][[j]],
    chaserPos[i][j][0] == chaserOrigPos[[i, j]]
    },
   {i, chaserNumber}, {j, 2}]]

Now the system is ready for NDSolve. Please note how WhenEvent stops the integration when any chaser is close enough ($\text{distance}<0.1$) to the target:

Manipulate[
 DynamicModule[{sol, tmax = 20},
  sol = With[{stopCondition = Or @@ (
        Norm[
            targetPos[targetVelocity, θ, t] - Through[(chaserPos[#] /@ {1, 2})[t]]
            ] < .1 & /@
         Range[chaserNumber])},

    NDSolve[{eqsys[chaserVelocity, targetVelocity, θ],
       WhenEvent[stopCondition, tmax = t; "StopIntegration"]
       } // Flatten,
     Flatten[
      Outer[chaserPos[#1][#2][#3] &, 
       Range[chaserNumber], {1, 2}, {t}]],
     {t, 0, tmax}]
    ];

  ParametricPlot[Evaluate[Append[
       Through[(chaserPos[#] /@ {1, 2})[t]] & /@ 
         Range[chaserNumber] /. sol[[1]],
       human[targetVelocity, θ, t]
       ]], {t, 0, tmax},
     Frame -> True, Axes -> False, ImageSize -> 400,
     Prolog -> {GrayLevel[.9], Polygon[chaserOrigPos]}
     ] /. Line[pts_] :> Arrow[pts] // Column[{tmax, #}] &
  ],

 {{θ, π/chaserNumber}, 0, (2 π)/chaserNumber},
 {{chaserVelocity, 10}, 0, 10},
 {{targetVelocity, 6}, 0, 10},
 ContinuousAction -> False]

chasing trail

And we can plot the survival time for different escaping directions:

With[{chaserVelocity = 10, targetVelocity = 9.5},
 tmaxLst = {};
  Table[
  Module[{tmax = 20},
   With[{stopCondition = Or @@ (
        Norm[
            targetPos[targetVelocity, θ, t] - Through[(chaserPos[#] /@ {1, 2})[t]]
            ] < .1 & /@
         Range[chaserNumber])},

    NDSolve[{eqsys[chaserVelocity, targetVelocity, θ],
       WhenEvent[stopCondition, tmax = t; "StopIntegration"]
       } // Flatten,
     Flatten[
      Outer[chaserPos[#1][#2][#3] &, 
       Range[chaserNumber], {1, 2}, {t}]],
     {t, 0, tmax}];

    AppendTo[tmaxLst, {θ, tmax}]
    ]],
  {θ, Rescale[Range[0, 1, .01], {0, 1}, {0, (2 π)/chaserNumber}]}];

 ListPolarPlot[tmaxLst,
  Joined -> True, PlotStyle -> Lighter[Red, .5],
  Mesh -> All, MeshStyle -> Red,
  Epilog -> {Gray,
    Line[{{0, 0}, Max[tmaxLst[[All, 2]]] {Cos[#], Sin[#]}}] & /@
     Rescale[Range[0, 1, 1/8], {0, 1}, {0, (2 π)/chaserNumber}],
    Circle[{0, 0}, #, {0, (2 π)/chaserNumber}] & /@
     FindDivisions[{0, Max[tmaxLst[[All, 2]]]}, 20]
    }]     ]

survival time curve

share|improve this answer

could anyone explain why my code fails? Namely, why does Mathematica think that raptor[[i]] is a scalar valued function even though its initial condition is clearly vector valued.

raptor[[i]] is a scalar valued element of something which you wish to be a vector valued object. You haven't told Mathematica that it is a vector. Whenever you call something which you wish to index in this way you must initialise it first.

As @Silvia mentioned in the comments, you should initialise raptor to be a list. Then you can call its elements by index i.

share|improve this answer

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.