How to solve a system of three coupled vectorized second order ODEs

Hi,

I am trying calculate the path of three stars in 3 dimensions using the ode45 function.

The formula goes as follow:

;

;

;

K1 and m are constants, r1, r2, r3 and v1 v2 v3 are vectors of the form . This means that I have three ode2 with six initial conditions in vector form. My latest attemp looks like this:

ode1 = @(t,r1,r2,r3) K1*(M2*(r2-r1)/norm(r2-r1)^3 + M3*(r3-r1)/norm(r3-r1)^3);

ode2 = @(t,r1,r2,r3) K1*(M1*(r1-r2)/norm(r2-r1)^3 + M3*(r3-r2)/norm(r3-r2)^3);

ode3 = @(r,r1,r2,r3) K1*(M1*(r1-r3)/norm(r3-r1)^3 + M2*(r2-r3)/norm(r3-r2)^3);

ode4 = @(t,v1,v2,v3) K2*v1;

ode5 = @(t,v1,v2,v3) K2*v2;

ode6 = @(t,v1,v2,v3) K2*v3;

odes = {ode1;ode2;ode3;ode4;ode5;ode6};

cond = [r1i,r2i,r3i,v1i,v2i,v3i];

span = linspace(0,20,500);

[tout, Sol] = ode45(odes,span,cond);

I tried using symbolic math as well but didn't get any further.

[Torsten]

How does r_1 look ? Is it r_1 = (r_11,r_12,r_13) (same for the other vectors) ?

[Samuel Gosselin]

Yes, here are my current r_i

r1i = [-0.5,0,0];

r2i = [0.5,0,0];

r3i = [0,1,0];

[Torsten]

But the r_i are the solution variables - how can you fix them ?

Or does it mean that acceleration is constant ?

Maybe for clarity you could write out the equations to be solved componentwise.

[Samuel Gosselin]

The r_i are position vectors for the different masses. I want to konw the trajectory described by the interaction of the masses. The three masses (m1, m2, m3) start at position r1i, r2i r3i with v1i, v2i, v3i velocities respectively.

I want to extract the values ri_11, ri_12, ri_13 to plot in a x,y,z coordinates system.

They do not have an initial acceleration but their interaction will produce an acceleration as the system evovle.

[David Goodmanson]

Hi Samuel,

You can make use of the fact that, since there are no external forces on the system, its center of mass moves with constant velocity. There is no need to keep that velocity around and it detracts from the solution. Without loss of generality you can set the center of mass to be at the origin, and set the momentum of the center of mass to be zero. Then

m1*r1 + m2*r2 + m3*r3 = 0    -->   r3 =  - (1/m3)*(m1*r1 + m2*r2)
m1*v1 + m2*v2 + m3*v3 = 0    -->   v3 =  - (1/m3)*(m1*v1 + m2*v2)

plugging those results into the equations eliminates r3,v3 and reduces the number of variables from 18 to 12. This probably works out best if m3 is the largest of the three masses.

It's also true that the total angular momentum of the system is conserved, but implementing that would be more difficult.

Answer 1

[Torsten]

System to solve for r1=(r1x,r1y,r1z) is:

dr1x/dt = v1x
dv1x/dt = K1*(m2*(r1x-r2x)/((r1x-r2x)^2+(r1y-r2y)^2+(r1z-r2z)^2)^1.5 + m3*(r1x-r3x)/((r1x-r3x)^2+(r1y-r3y)^2+(r1z-r3z)^2)^1.5)
dr1y/dt = v1y
dv1y/dt = K1*(m2*(r1y-r2y)/((r1x-r2x)^2+(r1y-r2y)^2+(r1z-r2z)^2)^1.5 + m3*(r1y-r3y)/((r1x-r3x)^2+(r1y-r3y)^2+(r1z-r3z)^2)^1.5)
dr1z/dt = v1z
dv1z/dt = K1*(m2*(r1z-r2z)/((r1x-r2x)^2+(r1y-r2y)^2+(r1z-r2z)^2)^1.5 + m3*(r1z-r3z)/((r1x-r3x)^2+(r1y-r3y)^2+(r1z-r3z)^2)^1.5)

Add the analogue equations for r2 and r3, give 18 initial conditions and use ode45 to solve.

A symbolic solution as you tried is impossible.