Electric field and a misbehaving particle

What I wish to do is to model the trajectory of a point particle of small negative charge in the field of a point charge of large positive mass, so the latter is stationary.

What I have done so far always gives the solution where the small particle is totally unaffected by the field it's in -- so it just travels on in a straight line, which is obviously rubbish. (It should be going in a circle centered at the origin!)

Could anyone please help me? Note that the code is in its rough/draft stage so signs and constants are not guaranteed, but that should not affect the general physics behind it!

We shall consider motion in a plane. The coordinates of the particle with small charge is (a(t),b(t)) and that with large charge is at (0,0).

function EField

clear 

trange = (0:.005:5);

%Starting conditions
a0 = ... ;
b0 = ... ;
adot0 = ... ;
bdot0 = ... ;
addot0 = -a0/((a0^2+b0^2)^(3/2));
bddot0 = -b0/((a0^2+b0^2)^(3/2));
vec0 = [a0, b0, adot0, bdot0, addot0, bddot0];

options = odeset('RelTol', 1e-8 ,'AbsTol', 1e-10);

[t,vec] = ode45(@EF, trange, vec0, options);

%Plot trajectory
plot(vec(:,1), vec(:,2));

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

This is the function EF in a separate m-file.

function sth = EF(t,vec)

a = vec(1);
b = vec(2);
adot = vec(3);
bdot = vec(4);
addot = vec(5);
bddot = vec(6);

sth = zeros(size(vec));
sth(1) = a;
sth(2) = b;
sth(3) = adot;
sth(4) = bdot;
sth(5) = -a/((a^2+b^2)^(3/2));
sth(6) = -b/(2*a^2+b^2)^(3/2));

Thank you.

[Jan]

What is "X" in "sth = zeros(size(X));"?!

[Richard]

Ah, that was a typo... For some reason I couldn't copy and paste, so I retyped the code -- unfortunately not too accurately!

Answer 1

[Jan]

What about:

...
sth(1) = adot;
sth(2) = bdot;
sth(3) = addot;
sth(4) = bddot;
...

The change of the position should be related to the velocity.

[Richard]

So there are only 4 entries for "sth"?

[Richard]

I may be misunderstood the mechanics of ode45. But I thought that vec and sth would specify the system of differential equations, so each entry in vec should correspond to an entry in sth. In a way it is something like vec is the LHS and sth is the RHS.

[Jan]

"sth" has 6 still components. The LHS of the function to be integrated contains the derivatives. So for the component which represents the position, the velocity is expected. And the velocities are changed by the accelerations. And the accelerations are changed by the external forces - in general.

[Richard]

Thank you Jan :) I am sure it is just me being silly, but I still don't quite understand what is going on. If

sth(1) = adot;

sth(2) = bdot;

sth(3) = addot;

sth(4) = bddot;

Then what are

sth(5) sth(6)?

And doesn't sth(1) = adot is equivalent to saying that LHS: vec(1)=a is equal to RHS sth(1)=adot? Or have I gotten the way ode45 works wrong?

[Jan]

I do not know to function you want to integrate. I assume, that sth(5) and (6) are correct in your original question and only the frist 4 components are wrong.

"sth" is the LHS, "vec" is the RHS of your function.

The function to be integrated is: dx/dt=... . It contains the local derivative of x. Let's use a simpler example: A point of mass 1 is accelerated in the gravitational field:

x''=g

The current state of the object is represented by x=[position, velocity]. Then the function to be integrated is:

function dx = dxdt(t, x)

dx(1) = x(2);

dx(2) = g;

[Richard]

Thanks so much, Jan. I have indeed misunderstood how ode45 works.