Advanced Search

Question: Problems with Solar System Model

Question:Problems with Solar System Model

nrebman1 15 Maple
dsolve numeric differential-equations
April 27 2016
1

Howdy all,

I am trying to create a solar system model by defining a force equation then using the sequence function to create a differential equation and then solving those differential equations using the initial conditions (in X,Y,Z coordinates). So far I have the code below.


> 

m[1] = 1.989*10^30; m[2] = 3.301*10^23; m[3] = 4.867*10^24; m[4] = 5.972*10^24+7.346*10^22; m[5] = 6.417*10^23; m[6] = 1.899*10^27; m[7] = 5.685*10^26; m[8] = 8.682*10^25; m[9] = 1.024*10^26; m[10] = 1.471*10^22; m[11] = 9.3*10^20; m[12] = 2.6*10^20; m[13] = 2*10^20; m[14] = 8.67*10^19; m[15] = 3.9*10^19; mass := Matrix(15, 1, [1.989*10^30, 3.301*10^23, 4.867*10^24, 5.972*10^24+7.346*10^22, 6.417*10^23, 1.899*10^27, 5.685*10^26, 8.682*10^25, 1.024*10^26, 1.471*10^22, 9.3*10^20, 2.6*10^20, 2*10^20, 8.67*10^19, 3.9*10^19]); G := 6.67408*10^(-11)

mass := Vector(4, {(1) = ` 15 x 1 `*Matrix, (2) = `Data Type: `*anything, (3) = `Storage: `*rectangular, (4) = `Order: `*Fortran_order})

  

0.6674080000e-10

 (1)
> 

sqrt(sum(x[i](t)^2, i = 1 .. 3));

(x[1](t)^2+x[2](t)^2+x[3](t)^2)^(1/2)

  

proc (i, j) options operator, arrow; sqrt(sum((x[i, k](t)-x[j, k](t))^2, k = 1 .. 3)) end proc

  

((x[1, 1](t)-x[3, 1](t))^2+(x[1, 2](t)-x[3, 2](t))^2+(x[1, 3](t)-x[3, 3](t))^2)^(1/2)

  

proc (i, j) options operator, arrow; [x[j, 1](t)-x[i, 1](t), x[j, 2](t)-x[i, 2](t), x[j, 3](t)-x[i, 3](t)] end proc

  

[x[3, 1](t)-x[1, 1](t), x[3, 2](t)-x[1, 2](t), x[3, 3](t)-x[1, 3](t)]

  

x[3, 1](t)-x[1, 1](t)

  

x[j, 1](t)-x[i, 1](t), x[j, 2](t)-x[i, 2](t), x[j, 3](t)-x[i, 3](t)

  

proc (i, j) options operator, arrow; [seq(x[j, k](t)-x[i, k](t), k = 1 .. 3)] end proc

  

[x[3, 1](t)-x[1, 1](t), x[3, 2](t)-x[1, 2](t), x[3, 3](t)-x[1, 3](t)]

  

x[3, 1](t)-x[1, 1](t)

 (2)
> 

diff(x[0, 1](t), t, t) = force(0)[1]:
> 

initialPositions := Matrix([[0, 0, 0], [-0.210e8, 0.426e8, 0.541e7], [0.106e9, -0.244e8, -0.644e7], [-0.139e9, -0.569e8, 0.316e4], [-0.177e9, -0.155e9, 0.111e7], [-0.802e9, 0.131e9, 0.174e8], [-0.480e9, -0.142e10, 0.438e8], [0.280e10, 0.103e10, -0.324e8], [0.420e10, -0.157e10, -0.645e8], [0.132e10, -0.477e10, 0.127e9], [0.431e9, -0.690e8, -0.816e8], [0.228e9, 0.305e9, -0.368e8], [0.300e9, -0.351e9, 0.217e9], [-0.434e9, -0.841e7, -0.284e8], [-0.115e9, -0.466e9, -0.612e8]])

initialPositions := Vector(4, {(1) = ` 15 x 3 `*Matrix, (2) = `Data Type: `*anything, (3) = `Storage: `*rectangular, (4) = `Order: `*Fortran_order})

 (3)
> 

initialVelocities := Matrix([[0, 0, 0], [-0.462e7, -0.170e7, 0.285e6], [0.666e6, 0.293e7, 0.183e4], [0.936e6, -0.239e7, 83.3], [0.145e7, -0.140e7, -0.650e5], [-0.196e6, -0.106e7, 0.879e4], [0.745e6, -0.271e6, -0.250e5], [-0.208e6, 0.524e6, 0.467e4], [0.161e6, 0.442e6, -0.129e5], [0.463e6, 0.294e5, -0.137e6], [0.193e6, 0.143e7, 0.923e4], [-0.119e7, 0.974e6, 0.116e6], [0.978e6, 0.562e6, -0.470e6], [0.166e6, -0.156e7, -0.131e5], [0.132e7, -0.170e6, 0.395e6]])

initialVelocities := Vector(4, {(1) = ` 15 x 3 `*Matrix, (2) = `Data Type: `*anything, (3) = `Storage: `*rectangular, (4) = `Order: `*Fortran_order})

 (4)
> 

ic1 := seq(seq(x[i, k](0) = initialPositions[i+1, k], k = 1 .. 3), i = 0 .. N-1); ic2 := seq(seq((D(x[i, k]))(0) = initialVelocities[i+1, k], k = 1 .. 3), i = 0 .. N-1); equations := {ic1, ic2, ode}; sol := dsolve(equations, numeric)

Error, (in dsolve/numeric/process_input) invalid specification of initial conditions, got {x[0, 1](0) = 0, x[0, 2](0) = 0, x[0, 3](0) = 0, x[1, 1](0) = -0.210e8, x[1, 2](0) = 0.426e8, x[1, 3](0) = 0.541e7, x[2, 1](0) = 0.106e9, x[2, 2](0) = -0.244e8, x[2, 3](0) = -0.644e7, x[3, 1](0) = -0.139e9, x[3, 2](0) = -0.569e8, x[3, 3](0) = 0.316e4, x[4, 1](0) = -0.177e9, x[4, 2](0) = -0.155e9, x[4, 3](0) = 0.111e7, x[5, 1](0) = -0.802e9, x[5, 2](0) = 0.131e9, x[5, 3](0) = 0.174e8, x[6, 1](0) = -0.480e9, x[6, 2](0) = -0.142e10, x[6, 3](0) = 0.438e8, x[7, 1](0) = 0.280e10, x[7, 2](0) = 0.103e10, x[7, 3](0) = -0.324e8, x[8, 1](0) = 0.420e10, x[8, 2](0) ...

  
> 

plots[odeplot](sol, [x[1, 1](t), x[1, 2](t), x[1, 3](t)], t = 0 .. 20, numpoints = 1000, axes = normal)

Error, (in plots/odeplot) input is not a valid dsolve/numeric solution

  
> 

plots[odeplot](sol, [seq(x[1, k](t), k = 1 .. 3)], t = 0 .. 20, numpoints = 1000, axes = normal)

Error, (in plots/odeplot) input is not a valid dsolve/numeric solution

  
> 

plots[odeplot](sol, [seq([seq(x[i, k](t), k = 1 .. 3)], i = 0 .. N-1)], t = 0 .. 20, numpoints = 1000, axes = normal)

Error, (in plots/odeplot) input is not a valid dsolve/numeric solution

  
> 

``


Download ass4.mw

Everything works fine until I try to execute the last line. When I do that I get an error that says "Error, in dsolve/numeric/process_input. invalid specifications of initial conditions.

At this point I am not sure if the problem lies in how I have defined my initial conditions or the way i've defined the force equation but I am open to any suggestions or ideas on where I should go from here.

Thanks in advance!

 

Nick
Please Wait...