# Euler's Method for 4-D Differential Equation (Computing)

1. Oct 15, 2007

### qspeechc

1. The problem statement, all variables and given/known data
Four dogs start running from the corners of a field 100m square. Label the dogs A, B, C, D clockwise in that order. Dog A chases dog B which chases dog C which chases dog D which chases dog A. Here 'chase' means that the direction in which the chasing dog moves is directly towards the dog it chases. All the dogs move at 7m/s. The purpose of this project is to compute their positions as a function of time, display paths, trajectories, and perhaps other interesting aspects.

2. Relevant equations
?

3. The attempt at a solution
This is from a MATLAB computing course I am taking. First of all, I am required to get the Differential equations (DE's) of the positions of the dogs. I have no idea how to do this, but I am thinking

d(pA)/dt = pB etc.
Forgive me, I do not have the time for latex. Here pX donates the position of dog X. Is this correct? ZHow do I break down the x- and y-components of position to give the x- and y- components of the velocity?

2. Oct 15, 2007

### D H

Staff Emeritus
That is not correct. For example say dog A starts at the upper right corner, dog B at the lower right corner. Dog A's initial position is (50,50) and dog B. (50,-50). You need dog A running toward dog B, not diagonally off the field.

You have been given the magnitude of each dog's velocity vector. You need to come up with a direction, and the problem statement tells you how to do that:

3. Oct 15, 2007

### qspeechc

Ok, if i want to break this down into x- and y- components of the position and velocity, it gets really hairy:

Update the x-component of the position of dog A by the previous x-component + the time step times the x-component of its velocity. And similarly for y-component and other dogs.

Update the x-component of the velocity of dog A by the previous x-component + the time step times ( B's x-component of position - A's x-component of position)/magnitude of the distance between the dogs, multiplied by 7. etc. for the rest.

Whew! what a calculation! And I need to program that junk! Is it correct?

4. Oct 15, 2007

### D H

Staff Emeritus
This is actually a very simple set of calculations. You only have four point objects and very limited interactions between the objects. Molecular simulations involve thousands of objects with many interactions, and forces as positions and velocities. Vehicle simulations add orientation to the whole mess.

5. Oct 15, 2007

### qspeechc

Ok, thanks but you're not helping. It is a 8-dimensional problem actually, two components for each dogs' position and velocity !

6. Oct 15, 2007

### D H

Staff Emeritus
So use vectors and matrices. Use loops. This is just a dozen (or maybe two dozen) lines of code.

7. Oct 15, 2007

### qspeechc

It's pretty hard to do in MATLAB, in VPython its a cinch. Anyway, I STILL do not know the way to do it using Euler's method, or the DE's for that matter. Would it be like the motion of two planets under the action of gravity, except they're not accelerating, and there are four dogs ( as opposed to two)?

8. Oct 15, 2007

### D H

Staff Emeritus
What makes this a cinch in VPython and not in Matlab? How could you solve this problem as specified in any language without building the differential equations?

9. Oct 16, 2007

### qspeechc

Does 90 lines of code seem appropriate for this problem? Those 90 lines STILL don't work. I am truely at my wit's end!

10. Oct 16, 2007

### robphy

If (as you say) it's easy in VPython, I suspect that you are taking advantage of the vector-algebraic functions like dot and norm, etc... I'm not sure if MATLAB has those features built in, but it shouldn't be hard to create those convenience functions.

Note that the dogs are accelerating... their velocity vectors are changing [although their magnitudes are constant].

11. Oct 16, 2007

### qspeechc

Yes, I meant it's easier to deal with vectors in VPython and all that.
Well, I have trimmed down my code, but my dogs all run off to the right! I am quite sure the initial conditions are correct.
Here is the code:

INITIALISATION OF PARAMETERS
a = sqrt((Bposx0-Aposx0)^2 + (Bposy0-Aposy0)^2);
aposx = Aposx0 + h*v*cos((Bposx0-Aposx0)/a);
Aposx = [Aposx aposx];
aposy = Aposy0 + h*v*sin(((Bposy0)-Aposy0)/a);
Aposy = [Aposy aposy];
ETC. SO ON FOR B AND C AND D

hold on

while ~(Aposx(i) == Bposx(i))&~(Aposy(i) == Bposy(i));
a = sqrt((Bposx(i)-Aposx(i))^2 + (Bposy(i)-Aposy(i))^2);
aposx = Aposx(i) + h*v*cos((Bposx(i)-Aposx(i))/a);
Aposx = [Aposx aposx];
aposy = Aposy(i) + h*v*sin((Bposy(i)-Aposy(i))/a);
Aposy = [Aposy aposy];
AND SIMILARLY FOR THE OTHER DOGGIES

if i>1
plot([Aposx(i-1) Aposx(i)], [Aposy(i-1) Aposy(i)], 'b')
plot([Bposx(i-1) Bposx(i)], [Bposy(i-1) Bposy(i)], 'g')
plot([Cposx(i-1) Cposx(i)], [Cposy(i-1) Cposy(i)], 'r')
plot([Dposx(i-1) Dposx(i)], [Dposy(i-1) Dposy(i)], 'k')
end
i = i+1;
end

12. Oct 16, 2007

### robphy

Why cos((Bposx0-Aposx0)/a) ? and similarly for sin?

13. Oct 16, 2007

### qspeechc

Because if we multiply the velocity with the cosine of the angle, we get the x-component, and the y-component if we multiply by the sine.

14. Oct 16, 2007

### robphy

Are you saying that "(Bposx0-Aposx0)/a" is the angle?

15. Oct 16, 2007

### qspeechc

Er, sorry I am not following you?

16. Oct 16, 2007

### D H

Staff Emeritus
Matrices and vectors are what Matlab excels at. You are not taking advantage of the features in Matlab. Since you are struggling, a bit of help here.

The obvious coordinate system for this problem has origin at the center of the field and axes aligned with the edges of the field. Suppose A is a 2-vector that represents the position of dog A and B is a 2-vector that represents the position of dog B. Placing dog A at the upper right corner and B at the lower right corner, the initial states of dogs A and B are
Code (Text):

A = [50 50];
B = [50 -50];

Dog A's velocity vector is in the direction of dog B, which means you need the unit vector from dog A to dog B:
Code (Text):

AtoB = B - A;
uhat = AtoB / sqrt(dot(AtoB, AtoB));

No trig functions at all here!

Since you have four dogs, it makes a bit more sense to represent the state (the positions of the dogs) as a 4x2 matrix rather than using four separate 2-vectors.

17. Oct 16, 2007

### qspeechc

Hey! No fair! I just figured that out on my own, right before you posted that :lol:!

Except your way is far more succinct, my code looked like this:

a = sqrt((Bpos(1)-Apos(1))^2 + (Bpos(2)-Apos(2))^2);
ra = [(Bpos(1)-Apos(1)) (Bpos(2)-Apos(2))];
Apos = Apos + h*v*ra./a;

Yuck! Yours looks far better! I never knew MATLAb could do dot products and so on. It really is extremely powerful, as I find out more and more every day.

By the way, I have this while loop, I want it to stop when the dog's are in the vicinity of one another, because they never actually touch (trust me, I've been doing it for time step of 1e-3, and have waited many minutes), is there any way I can do this? Currently my code is:

while ~(Apos == Bpos)
a = sqrt((Bpos(1)-Apos(1))^2 + (Bpos(2)-Apos(2))^2);

So it never actually stops. How do I get it to stop when they are near one another?

18. Oct 16, 2007

### robphy

So, you got rid of those spurious trig functions. Great.

In an iterative calcuation, it's not likely that Apos == Bpos will ever be achieved.
More practically, you can stop when the magnitude of the displacement (or better, the square-distance) is sufficiently small...
i.e. dot(AtoB,AtoB)<1.0e-8.

19. Oct 16, 2007

### D H

Staff Emeritus
robphy gave a very good answer here. In numerical computing, it is usually unwise to use == when looking for agreement between two floating point numbers. Remember this. A lot of engineers and scientists who do numerical computing don't remember this and end up having problems of some sort. Only use == if you really, really mean it. The best thing to do in most cases is to compare the square error. (That way you aren't taking a superfluous square root).

BTW, the dogs never do touch because each is following the same well-known mathematical curve, rotated by 90 degrees for each dog.

20. Oct 16, 2007

### qspeechc

Thanks! EDIT: what's a floating point number? We haven't gotten into that, and I do not think we will.

And, last question, I promise: can I, and how do you, solve this with MATLAB's built in ODE solvers? Can it even handle 8 dimensions?
I tried to define the function file:

function f= doggy(t, x)

f= zeros(4, 2);
f(1) = 7*(x(2) - x(1))/(sqrt(dot((x(2)-x(1)),(x(2)-x(1)))));
f(2) = 7*(x(3) - x(2))/(sqrt(dot((x(3)-x(2)),(x(3)-x(2)))));
f(3) = 7*(x(4) - x(3))/(sqrt(dot((x(4)-x(3)),(x(4)-x(3)))));
f(4) = 7*(x(1) - x(4))/(sqrt(dot((x(1)-x(4)),(x(1)-x(4)))));

the x(a)'s supposed to be the position vectors of the dogs. In the command window I entered:

x0 = [0 100 100 0; 100 100 0 0];
[t, x] = ode45(@doggy, [0 10], x0); plot(t,x)

x0 is supposed to give their initial positions, as I had defined them.

I get the incomprehensible (what else) error:

Warning: Divide by zero.
> In doggy at 5
In funfun\private\odearguments at 110
In ode45 at 173
Warning: Divide by zero.
> In doggy at 6
In funfun\private\odearguments at 110
In ode45 at 173
??? Error using ==> funfun\private\odearguments
DOGGY must return a column vector.

Error in ==> ode45 at 173
[neq, tspan, ntspan, next, t0, tfinal, tdir, y0, f0, odeArgs, ...

HEEEEELP!
EDIT AGAIN: how do you solve this analytically? Vector differential equation?

Last edited: Oct 16, 2007