Hi,
the following equations will produce some output.
sols = NDSolve[{x1''[
t] == (G m2 (x2[t] -
x1[t]))/((x2[t] - x1[t])^2 - (y2[t] - y1[t])^2)^(3/
2) + (G m3 (x3[t] -
x1[t]))/((x3[t] - x1[t])^2 - (y3[t] - y1[t])^2)^(3/2),
y1''[t] == (G m2 (y2[t] -
y1[t]))/((x2[t] - x1[t])^2 - (y2[t] - y1[t])^2)^(3/
2) + (G m3 (y3[t] -
y1[t]))/((x3[t] - x1[t])^2 - (y3[t] - y1[t])^2)^(3/2),
x2''[t] == (G m1 (x1[t] -
x2[t]))/((x1[t] - x2[t])^2 - (y1[t] - y2[t])^2)^(3/
2) + (G m3 (x3[t] -
x2[t]))/((x3[t] - x2[t])^2 - (y3[t] - y2[t])^2)^(3/2),
y2''[t] == (G m1 (y1[t] -
y2[t]))/((x1[t] - x2[t])^2 - (y1[t] - y2[t])^2)^(3/
2) + (G m3 (y3[t] -
y2[t]))/((x3[t] - x2[t])^2 - (y3[t] - y2[t])^2)^(3/2),
x3''[t] == (G m1 (x1[t] -
x3[t]))/((x1[t] - x3[t])^2 - (y1[t] - y3[t])^2)^(3/
2) + (G m2 (x2[t] -
x3[t]))/((x2[t] - x3[t])^2 - (y2[t] - y3[t])^2)^(3/2),
y3''[t] == (G m1 (y1[t] -
y3[t]))/((x1[t] - x3[t])^2 - (y1[t] - y3[t])^2)^(3/
2) + (G m2 (y2[t] -
y3[t]))/((x2[t] - x3[t])^2 - (y2[t] - y3[t])^2)^(3/2),
x1[0] == d1, y1[0] == 0, x1'[0] == -0.8, y1'[0] == v1,
x2[0] == d2, y2[0] == 0, x2'[0] == -0.2, y2'[0] == v2,
x3[0] == d3, y3[0] == 0, x3'[0] == -0.1,
y3'[0] == v3} /. {G -> 0.1, m1 -> 1, m2 -> 1, m3 -> 1, d1 -> 10,
v1 -> 1, d2 -> 20, v2 -> 2, d3 -> 30, v3 -> 1}, {x1[t], y1[ t],
x2[t], y2[t], x3[t], y3[t]}, {t, 0, 10}];
You can plot functions like this:
Plot[{x1[t], x2[t]} /. sols, {t, 0, 10}]
You can also use ParametricPlot:
ParametricPlot[{{x1[t], y1[t]} /. sols, {x2[t], y2[t]} /. sols, {x3[t], y3[t]} /. sols}, {t, 0, 10}, PlotRange -> {{-30, 30}, {-30, 30}}]
The results are still not what you would expect. This is because your equations for the gravitational force are not correct. Where you divide by the squared distances in the x and y direction you should actually add not subtract:
sols = NDSolve[{x1''[
t] == (G m2 (x2[t] -
x1[t]))/((x2[t] - x1[t])^2 + (y2[t] - y1[t])^2)^(3/
2) + (G m3 (x3[t] -
x1[t]))/((x3[t] - x1[t])^2 + (y3[t] - y1[t])^2)^(3/2),
y1''[t] == (G m2 (y2[t] -
y1[t]))/((x2[t] - x1[t])^2 + (y2[t] - y1[t])^2)^(3/
2) + (G m3 (y3[t] -
y1[t]))/((x3[t] - x1[t])^2 + (y3[t] - y1[t])^2)^(3/2),
x2''[t] == (G m1 (x1[t] -
x2[t]))/((x1[t] - x2[t])^2 + (y1[t] - y2[t])^2)^(3/
2) + (G m3 (x3[t] -
x2[t]))/((x3[t] - x2[t])^2 + (y3[t] - y2[t])^2)^(3/2),
y2''[t] == (G m1 (y1[t] -
y2[t]))/((x1[t] - x2[t])^2 + (y1[t] - y2[t])^2)^(3/
2) + (G m3 (y3[t] -
y2[t]))/((x3[t] - x2[t])^2 + (y3[t] - y2[t])^2)^(3/2),
x3''[t] == (G m1 (x1[t] -
x3[t]))/((x1[t] - x3[t])^2 + (y1[t] - y3[t])^2)^(3/
2) + (G m2 (x2[t] -
x3[t]))/((x2[t] - x3[t])^2 + (y2[t] - y3[t])^2)^(3/2),
y3''[t] == (G m1 (y1[t] -
y3[t]))/((x1[t] - x3[t])^2 + (y1[t] - y3[t])^2)^(3/
2) + (G m2 (y2[t] -
y3[t]))/((x2[t] - x3[t])^2 + (y2[t] - y3[t])^2)^(3/2),
x1[0] == d1, y1[0] == 0, x1'[0] == -0.8, y1'[0] == v1,
x2[0] == d2, y2[0] == 0, x2'[0] == -0.2, y2'[0] == v2,
x3[0] == d3, y3[0] == 0, x3'[0] == -0.1,
y3'[0] == v3} /. {G -> 2, m1 -> 1, m2 -> 1, m3 -> 1, d1 -> 1,
v1 -> 0.1, d2 -> 2, v2 -> 0.2, d3 -> 3, v3 -> 0.4}, {x1[t],
y1[ t], x2[t], y2[t], x3[t], y3[t]}, {t, 0, 100}];
You will notice that this is probably still not what you want, or is it?
Cheers,
Marco
