# Why this system of equations doesn't work [closed]

This system of equations is interesting because it should describe the motion of three bodies in the space. System doesn't work because of mistakes in code, I guess. Another problem is how to set conditions to satisfy periodic or almost periodic solution. I hope that $Mathematica$ has some functions to set the parameters on right way and to get periodic solution.

sol1 = NDSolveValue[{x1''[t] == (x2[t] - x1[t])/r12[t]^3 + (
x3[t] - x1[t])/r13[t]^3,
y1''[t] == (y2[t] - y1[t])/r12[t]^3 + (y3[t] - y1[t])/r13[t]^3,
z1''[t] == (z2[t] - z1[t])/r12[t]^3 + (z3[t] - z1[t])/r13[t]^3,
x2''[t] == (x1[t] - x2[t])/r12[t]^3 + (x3[t] - x2[t])/r23[t]^3,
y2''[t] == (y1[t] - y2[t])/r12[t]^3 + (y3[t] - y2[t])/r23[t]^3,
z2''[t] == (z1[t] - z2[t])/r12[t]^3 + (z3[t] - z2[t])/r23[t]^3,
x3''[t] == (x1[t] - x3[t])/r13[t]^3 + (x2[t] - x3[t])/r23[t]^3,
y3''[t] == (y1[t] - y3[t])/r13[t]^3 + (y2[t] - y3[t])/r23[t]^3,
z3''[t] == (z1[t] - z3[t])/r13[t]^3 + (z2[t] - z3[t])/r23[t]^3,
y1[t]*z1'[t] - z1[t]*y1'[t] + y2[t]*z2'[t] - z2[t]*y2'[t] +
y3[t]*z3'[t] - z3[t]*y3'[t] = 0,
x1[t]*z1'[t] - z1[t]*x1'[t] + x2[t]*z2'[t] - z2[t]*x2'[t] +
x3[t]*z3'[t] - z3[t]*x3'[t] = 0,
x1[t]*y1'[t] - y1[t]*x1'[t] + x2[t]*y2'[t] - y2[t]*x2'[t] +
x3[t]*y3'[t] - y3[t]*x3'[t] = 0,
r12[t] =
Abs[((x2[t] - x1[t])^2 + (y2[t] - y1[t])^2 + (z2[t] - z1[t])^2)]^(
1/2),
r13[t] =
Abs[((x3[t] - x1[t])^2 + (y3[t] - y1[t])^2 + (z3[t] - z1[t])^2)]^(
1/2), r23[t] =
Abs[((x3[t] - x2[t])^2 + (y3[t] - y2[t])^2 + (z3[t] - z2[t])^2)]^(
1/2), x1[0] == 1, y1[0] == 1, z1[0] == 1, x1'[0] == 1,
y1'[0] == 1, z1'[0] == 1, x2[0] == 2, y2[0] == 2, z2[0] == 2,
x2'[0] == 2, y2'[0] == 2, z2'[0] == 2, x3[0] == 3, y3[0] == 3,
z3[0] == 3, x3'[0] == 3, y3'[0] == 3, z3'[0] == 3}, {x1, y1, z1,
x2, y2, z2, x3, y3, z3}, {t, 0, 1}]

-
Try naming NDSolve instead of NDSolveValue. Also there are three spots where you have used = (Set) instead of == (Equal). –  Thies Heidecke Mar 17 '13 at 14:15
yes, but, how to plot three functions in 3d and how to use == in case on right hand side have zero –  Pipe Mar 17 '13 at 14:33
try ParametricPlot3D[{x1[t], y1[t], z1[t]} /. sol1, {t, 0, 1}] –  Thies Heidecke Mar 17 '13 at 14:36
@ThiesHeidecke NDSolveValue is new in v9. –  rcollyer Mar 17 '13 at 15:13
@rcollyer oops! thanks for pointing that out, i should catch up on the new functions of v9 i guess... –  Thies Heidecke Mar 17 '13 at 15:42
show 1 more comment

## closed as too localized by Yves Klett, ssch, m_goldberg, Oleksandr R., halirutanMar 17 '13 at 23:21

This question is unlikely to help any future visitors; it is only relevant to a small geographic area, a specific moment in time, or an extraordinarily narrow situation that is not generally applicable to the worldwide audience of the internet. For help making this question more broadly applicable, visit the help center.If this question can be reworded to fit the rules in the help center, please edit the question.

There are some instances of = (Set) where you meant == (Equal) in your code. Also you are using three helper functions r12, r13, r23, that don't show up in the final solution. In that case you can incorporate them in your original equations by Replaceing (/.) their occurences with your definitions, or you could choose to let them be part of your output variables. Here's the first of the two options:

helper = {
r12[t] -> Abs[((x2[t] - x1[t])^2 + (y2[t] - y1[t])^2 + (z2[t] - z1[t])^2)]^(1/2),
r13[t] -> Abs[((x3[t] - x1[t])^2 + (y3[t] - y1[t])^2 + (z3[t] - z1[t])^2)]^(1/2),
r23[t] -> Abs[((x3[t] - x2[t])^2 + (y3[t] - y2[t])^2 + (z3[t] - z2[t])^2)]^(1/2)
};
odes = {
x1''[t] == (x2[t] - x1[t])/r12[t]^3 + (x3[t] - x1[t])/r13[t]^3,
y1''[t] == (y2[t] - y1[t])/r12[t]^3 + (y3[t] - y1[t])/r13[t]^3,
z1''[t] == (z2[t] - z1[t])/r12[t]^3 + (z3[t] - z1[t])/r13[t]^3,
x2''[t] == (x1[t] - x2[t])/r12[t]^3 + (x3[t] - x2[t])/r23[t]^3,
y2''[t] == (y1[t] - y2[t])/r12[t]^3 + (y3[t] - y2[t])/r23[t]^3,
z2''[t] == (z1[t] - z2[t])/r12[t]^3 + (z3[t] - z2[t])/r23[t]^3,
x3''[t] == (x1[t] - x3[t])/r13[t]^3 + (x2[t] - x3[t])/r23[t]^3,
y3''[t] == (y1[t] - y3[t])/r13[t]^3 + (y2[t] - y3[t])/r23[t]^3,
z3''[t] == (z1[t] - z3[t])/r13[t]^3 + (z2[t] - z3[t])/r23[t]^3
(* left out definitions below to avoid overdetermined system *)
(*,
y1[t]*z1'[t]-z1[t]*y1'[t]+y2[t]*z2'[t]-z2[t]*y2'[t]+y3[t]*z3'[t]-z3[t]*y3'[t]==0,
x1[t]*z1'[t]-z1[t]*x1'[t]+x2[t]*z2'[t]-z2[t]*x2'[t]+x3[t]*z3'[t]-z3[t]*x3'[t]==0,
x1[t]*y1'[t]-y1[t]*x1'[t]+x2[t]*y2'[t]-y2[t]*x2'[t]+x3[t]*y3'[t]-y3[t]*x3'[t]==0
*)
} /. helper
ics = {
x1[0]  == 1, y1[0]  == 1, z1[0]  == 1, x1'[0] == 1, y1'[0] == 1, z1'[0] == 1,
x2[0]  == 2, y2[0]  == 2, z2[0]  == 2, x2'[0] == 2, y2'[0] == 2, z2'[0] == 2,
x3[0]  == 3, y3[0]  == 3, z3[0]  == 3, x3'[0] == 3, y3'[0] == 3, z3'[0] == 3
}

sol1 = NDSolve[ Join[odes, ics], {x1, y1, z1, x2, y2, z2, x3, y3, z3}, {t, 0, 1}]


Also i left out three of your differential equations because the system is already fully determined with the first 9 equations.

-
thank you very much for this post..I just need to plod 3 functions in 3d ... hope that it will work.. do you have idea how to satisfy periodic solution? –  Pipe Mar 17 '13 at 14:41
You can plot multiple functions in 3D with ParametricPlot3D. Also have a look at Introduction to Advanced Numerical Differential Equation Solving in Mathematica. The documentation also has an example regarding periodic boundary conditions. –  Thies Heidecke Mar 17 '13 at 14:52
thank you very much Thies, helpful post and comments –  Pipe Mar 17 '13 at 15:12