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

[qspeechc]

Homework Statement

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.

Homework Equations

?

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?

 

[D H]

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?

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:

Here 'chase' means that the direction in which the chasing dog moves is directly towards the dog it chases.

 

[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?

 

[D H]

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.

 

[qspeechc]

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

 

[D H]

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

 

[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)?

 

[D H]

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?

 

[qspeechc]

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

 

[robphy]

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

post portions of your code

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)?

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].

 

[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

 

[robphy]

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

 

[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.

 

[robphy]

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

 

[qspeechc]

Er, sorry I am not following you?

 

[D H]

Yes, I meant it's easier to deal with vectors in VPython and all that.

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

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:

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.

 

[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?

 

[robphy]

Apos = Apos + h*v*ra./a;

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

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?

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.

 

[D H]

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 ...

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.

 

[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?

 

[D H]

To use the built-in ODE solvers you have to work with vectors. You are trying to use matrices (and making some errors in doing so). The typical approach is to pile the various state elements together. For example, state = [ dog A x pos, dog A y pos, dog B x pos, ... ]. The only limit is the size of memory. A lot of times different elements of the state "vector" will have different units. You have a simple case; the state is just position.

The 'reshape' command provides another option when the state is more naturally represented in matrix form, as is the case here. You form the initial state as a 4x2 matrix and then reshape it to 8x1 for compatibility with the built-in functions. In the derivative function, reshape the 8x1 state to 4x2, compute the derivatives as a 4x2 matrix, and then reshape to 8x1 for output.

I mentioned you had some errors. For example,

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

does not do what you want it to do. Suppose here f and x are each 4x2. What you meant to say was

f(1,:) = 7*(x(2,:) - x(1,:))/(sqrt(dot((x(2,:)-x(1,:)),(x(2,:)-x(1,:)))));

The colon is very important. Type 'help colon' in your Matlab command window.

 

[qspeechc]

Er, that completely went over my head! Probably because we haven't covered the built in ODE solvers yet. I am turning in and will continue with this tomorrow, hopefully.

Thank you very much D H and robphy! You guys (I assume you're male, sorry, slap on the wrists for me!) have been a great help. Thanks for all the time and insight you put into this, even though this problem is beneath you! Thank you very much!

 

[D H]

Do try using the command 'help colon'. Matlab's help files are quite helpful and usually well written. Guys is correct in my case. No need for the apology though, as 'guys' has become a unisex term.

 

[qspeechc]

Sorry to harp on about this!
Ok, I defined dogA to be at [0 100] initially, dogB at [100 100] etc.
So trying to solve this with ode45, I wrote the function file doggy.m
where x(1) is the x-position of dogA, x(2) is the y-position of dogA, x(3) the x-position of dogB etc.:

function f= doggy(t, x)

f= zeros(8, 1);
f(1) = 7*(x(3) - x(1))/(sqrt((x(3)-x(1))^2 + (x(4)-x(2))^2));
f(2) = 7*(x(4) - x(2))/(sqrt((x(3)-x(1))^2 + (x(4)-x(2))^2));
f(3) = 7*(x(5) - x(3))/(sqrt((x(5)-x(3))^2 + (x(6)-x(4))^2));
f(4) = 7*(x(6) - x(4))/(sqrt((x(5)-x(3))^2 + (x(6)-x(4))^2));
f(5) = 7*(x(7) - x(5))/(sqrt((x(7)-x(5))^2 + (x(8)-x(6))^2));
f(6) = 7*(x(8) - x(6))/(sqrt((x(7)-x(5))^2 + (x(8)-x(6))^2));
f(7) = 7*(x(1) - x(7))/(sqrt((x(1)-x(7))^2 + (x(2)-x(8))^2));
f(8) = 7*(x(2) - x(8))/(sqrt((x(1)-x(7))^2 + (x(2)-x(8))^2));

Then I made a separate script file:

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

But it gives a horribly wrong plot; two dogs start at [0 100] and two at [0 0], and they run off to the right. Is x0 wrong? I thought x0 as I have it gives the
[Ax(0) Ay(0) Bx(0) ... Dy(0)] corresponding to the dogs' initial positions.

 

[qspeechc]

I have no clue any more what's going on! I tried this instead, because the above code does not correct the velocity, it does not resolve it into the x- and y- components, so I tried this:

function f= doggy(t, x)

v= 7;

f(4, 2) = v*([x(1)/sqrt(dot(x(1),x(1))); x(2)/sqrt(dot(x(2),x(2)));...
x(3)/sqrt(dot(x(3),x(3))); x(4)/sqrt(dot(x(4),x(4)))])
f= zeros(reshape(f(4,2),[8,1]));
resphape(reshape(f(4,2),[8,1]))

And ofcourse I got an error:

? Error using ==> run
Input argument "x" is undefined.

And besides, if I got that code to work, I still don't know how to plot it!

I am starting to think this cannot be solved with MATLAB's built in ODE solvers, because it MUST be treated in matrix form.

 

[qspeechc]

Last effort:

function f= doggy(t, y)

v= 7;

f= zeros(8,1);
g= reshape(f, [2,4]);
x= reshape(y, [2,4]);

h1= (sqrt(dot((x(:,2) - x(:,1)),(x(:,2) - x(:,1)))));
g(:,1)= v*(x(:,2) - x(:,1))/h1;
h2= (sqrt(dot((x(:,3) - x(:,2)),(x(:,3) - x(:,2)))));
g(:,2)= v*(x(:,3) - x(:,2))/h2;
h3= (sqrt(dot(x(:,4) - x(:,3)),(x(:,4) - x(:,3)))));
g(:,3)= v*(x(:,4) - x(:,3))/h3;
h4= (sqrt(dot(x(:,1) - x(:,4)),(x(:,1) - x(:,4)))));
g(:,4)= v*(x(:,1) - x(:,4))/h4;

y= reshape(x, [8,1]);
f= reshape(g, [8,1]);

And ofcourse the error is:

? Error using ==> run
Input argument "y" is undefined.

How do I go about defining y? And then plotting this? Can't be the same script file as in post #24, unless the code I have written has a 8x1 output. HEEEEEEELP! I'm confused =(.

 

[qspeechc]

Ok^ that, I think, is fine, but plotting the results is another story:

x0 = [0 100 100 100 100 0 0 0];
[t,y]=ode45(@doggy,[0 30], x0);
plot(y(1),y(2),'b',y(3),y(4),'g',y(5),y(6),'r',y(7),y(8),'k')

Which is supposed to plot the positions and all that, but it just plots one point, very weird...