# How to solve a system of differential equations with no exact solution without specifying all of the initial conditions

I have a system of differential equations that does not have an exact solution to them. I get around this by using ndsolvevalue which works great if you provide actual numbers in the initial conditions. However, what I want to do now is leave some variables unspecified so that, upon solving the differential equations, the end result is an equation that is a function of these unspecified variables. I have included the actual code beneath. Basically, I'm trying to solve the system of diff eqs without giving a value to m1 or m2. Is this possible in Mathematica or must the initial conditions have a number attached to them (i.e. are not symbolic in nature)?

G = 6.6743*10^-11;
c = 2.99792458*10^8;
Periods = 5.599829;
P = Periods*86400;
ecc = 0.018;
a = Power[(P^2*G*(m1 + m2))/(4*π^2), (3)^-1]
μ = (m1 m2)/(m1 + m2);
aph = a (1 + ecc);
x0 = (m2*aph)/(m1 + m2)
y0 = 0;
X0 = -((m1*aph)/(m1 + m2))
Y0 = 0;
Vel1 = ((G μ )/a  m2/m1 (1 - ecc)/(1 + ecc) )^(
1/2)(*initial velocity star 1*)
Vel2 = -((G μ )/a  m2/m1 (1 - ecc)/(1 + ecc) )^(1/2) m1/
m2(*initial velocity star 2*)

system1 = {
x1''[t] == -((G m1 m2  (x1[t] - x2[t]))/(
m1 ((x1[t] - x2[t])^2 + (y1[t] - y2[t])^2)^(3/2))),
y1''[t] == -((G m1 m2  (y1[t] - y2[t]))/(
m1 ((x1[t] - x2[t])^2 + (y1[t] - y2[t])^2)^(3/2))),
x2''[t] == (G m1 m2  (x1[t] - x2[t]))/(
m2 ((x1[t] - x2[t])^2 + (y1[t] - y2[t])^2)^(3/2)),
y2''[t] == (G m1 m2  (y1[t] - y2[t]))/(
m2 ((x1[t] - x2[t])^2 + (y1[t] - y2[t])^2)^(3/2))};

initials = {
x1[0] == x0,
x1'[0] == 0,
y1[0] == 0,
y1'[0] == Vel1,
x2[0] == X0,
x2'[0] == 0,
y2[0] == 0,
y2'[0] == Vel2};

{xbh1, ybh1, xbh2, ybh2} =
DSolve[{system1, initials}, {x1, y1, x2, y2}, t];
bh1 = {(xbh1[t*86400]/(1.5*10^11)), (ybh1[
t*86400]/(1.5*10^11))};(*larger star position vectors*)
bh2 = {(xbh2[t*86400]/(1.5*10^11)), (ybh2[
t*86400]/(1.5*10^11))};(*smaller star position vectors*)

• Welcome to Mathematica SE. To get started:1) take the introductory tour now,2) when you see good questions and answers, vote them up by clicking the gray triangles, because the credibility of the system is based on the reputation gained by users sharing their knowledge,3) remember to accept the answer, if any, that solves your problem, by clicking checkmark sign,4) give help too, by answering questions in your areas of expertise.
– kcr
Mar 14 at 3:01
• At a minimum, use center-of-mass coordinates to eliminate two dependent variables. FullSimplify[SubtractSides[system1[[1]], system1[[3]]] /. {x2 -> Function[{t}, x1[t] + x[t]], y2 -> Function[{t}, y1[t] + y[t]]}], and similarly for the other two equations. Mar 14 at 23:03

I'm trying to solve the system of diff eqs without giving a value to m1 or m2. Is this possible in Mathematica

Use ParametricNDSolve like this

ParametricNDSolve[{system1, initials},{x1[t],y1[t],x2[t],y2[t]},{t,0,1},{m1,m2}]


the end result is an equation that is a function of these unspecified variables"

you can't get an "equation" from numerical solver. Another option other than using ParametricNDSolve is to use Manipulate. Make a slider for m1 and slider for m2. Each time you change the slider, you call NDSolve with the current values for m1 and m2.

• (+1) and sorry. I only saw after I posted
– kcr
Mar 14 at 2:55
• Thank you for the help! I was not aware of ParametricNDSolve until now! Yes, equation was a poor choice of words on my part. There is one issue I'm still finding however. These diff eqs will be used to fit a model to some data. This was not an issue when I specified m1 and m2 (my nonlinearmodelfit worked great). However, now that m1 and m2 are a part of these Parametric Functions and are treated as fitting parameters, my Nonlinearmodel fit appears to be having a fit! I get an error that the function value is not a list of real numbers. Any thoughts? Mar 14 at 4:33
• "get an error that the function value is not a list of real numbers. Any thoughts?" may be you are not using the solution from ParametricNDSolve correctly in Nonlinearmodel fit. It would be better to post separate question on this specific issue with all the code you used and the error you got. Mar 14 at 8:19
• Good idea! Thank you again for the help! Mar 14 at 13:50
• Hello again! If my previous question about nonlinearmodelfit intrigued you, I have made it a separate post here: mathematica.stackexchange.com/questions/265112/… Mar 14 at 18:39

G = 6.6743*10^-11;
c = 2.99792458*10^8;
Periods = 5.599829;
P = Periods*86400;
ecc = 0.018;
a = Power[(P^2*G*(m1 + m2))/(4*\[Pi]^2), (3)^-1]
\[Mu] = (m1 m2)/(m1 + m2);
aph = a (1 + ecc);
x0 = (m2*aph)/(m1 + m2)
y0 = 0;
X0 = -((m1*aph)/(m1 + m2))
Y0 = 0;
Vel1 = ((G \[Mu])/a m2/m1 (1 - ecc)/(1 + ecc))^(1/
2)(*initial velocity star 1*)
Vel2 = -((G \[Mu])/a m2/
m1 (1 - ecc)/(1 + ecc))^(1/2) m1/m2(*initial velocity star 2*)

system1 = {x1''[
t] == -((G m1 m2 (x1[t] -
x2[t]))/(m1 ((x1[t] - x2[t])^2 + (y1[t] - y2[t])^2)^(3/
2))), y1''[
t] == -((G m1 m2 (y1[t] -
y2[t]))/(m1 ((x1[t] - x2[t])^2 + (y1[t] - y2[t])^2)^(3/
2))), x2''[
t] == (G m1 m2 (x1[t] -
x2[t]))/(m2 ((x1[t] - x2[t])^2 + (y1[t] - y2[t])^2)^(3/2)),
y2''[t] == (G m1 m2 (y1[t] -
y2[t]))/(m2 ((x1[t] - x2[t])^2 + (y1[t] - y2[t])^2)^(3/2))};

initials = {x1[0] == x0, x1'[0] == 0, y1[0] == 0, y1'[0] == Vel1,
x2[0] == X0, x2'[0] == 0, y2[0] == 0, y2'[0] == Vel2};


and then we need to use ParametricNDSolve to have some unspecified parameters.

ParametricNDSolve[{system1, initials}, {x1[t], y1[t], x2[t],
y2[t]}, {t, 0, 1}, {m1, m2}]


which runs without any issues