Message Boards Message Boards

0
|
3840 Views
|
2 Replies
|
0 Total Likes
View groups...
Share
Share this post:
GROUPS:

Interpolating Functions

Posted 10 years ago

Hey everyone,

I am currently working on a senior project at my university and the project is about using my code to give a rough estimate where the Pioneer 10 spacecraft is. I used Newton's force equations and set up a numerical three-body diagram problem with NDSolve to get a position for the spacecraft. I notice that when I run this I get an interpolating function with a certain value. Am I then supposed to but that value in for the function relating to the spacecraft for a position?

Thank you

Brandyn Wessels

POSTED BY: Brandyn Wessels
2 Replies
Posted 10 years ago

Hello David,

This is my code for the trajectory of the Pioneer 10 spacecraft using a small solar system, (Sun, Earth, and Jupiter) and only accounting for gravitational forces. I'm working in AU units and trying to use the correct initial conditions for position and velocity at launch to the location of the craft in 1987. That's why the time of integration stops at 15 years after launch. So for the project I want to make sure I get accurate coordinates in the (x,y,z) plane from the Sun and compare my calculation with measured data.

In[1]:=
Clear[r1,r2,r3,r4,m1,m2,m3,m4,G,M,x,y,z,t]
r1={x1[t],y1[t],z1[t]};
r2={x2[t],y2[t],z2[t]};
r3={x3[t],y3[t],z3[t]};
r4={x4[t],y4[t],z4[t]};
r12=r2-r1;
r13=r3-r1;
r14=r4-r1;
r23=r3-r2;
r21=r1-r2;
r24=r4-r2;
r31=r1-r3;
r32=r2-r3;
r34=r4-r3;
r41=r1-r4;
r42=r2-r4;
r43=r3-r4;

equation1= D[r1,{t,2}]-  G m2 r12/(r12.r12)^(3/2+\[Delta])-
     m3 G r13/(r13.r13)^(3/2+\[Delta])- G m4 r14/(r14.r14)^(3/2+\[Delta]);

equation2= D[r2,{t,2}]-m1 G r21/(r21.r21)^(3/2+\[Delta])-
    m3 G r23/(r23.r23)^(3/2+\[Delta])- G m4 r24/(r24.r24)^(3/2+\[Delta]);

equation3= D[r3,{t,2}]- m2 G r32/(r32.r32)^(3/2+\[Delta])-
    m1 G r31/(r31.r31)^(3/2+\[Delta])- G m4 r34/(r34.r34)^(3/2+\[Delta]);

equation4= D[r4,{t,2}]-  G m2 r42/(r42.r42)^(3/2+\[Delta])-
     m3 G r43/(r43.r43)^(3/2+\[Delta])-    G m1 r41/(r41.r41)^(3/2+\[Delta]);

In[26]:=
Clear[M,m1,m2,m3,m4,G]
list1=Thread[(r1/.t->0)-{0,0,0}==0];
list2=Thread[(D[r1,t]/.t->0)-{0,0,0}==0];
list3=Thread[(r2/.t->0)-{1,0,0}==0];
list4=Thread[(D[r2,t]/.t->0)-{0,Sqrt[G M/x2[0]],0}==0];
list5=Thread[(r3/.t->0)-{4.95,0,0}==0];
list6=Thread[(D[r3,t]/.t->0)-{0,-Sqrt[G M/x3[0]],0}==0];
list7=Thread[(r4/.t->0)-{0.999940891,1.000041014,1.000000001}==0];
list8=Thread[(D[r4,t]/.t->0)-{0,(2.4512*10^-7),0}==0];
EquationList=Join[Thread[equation1==0],
    Thread[equation2==0],
    Thread[equation3==0],Thread[equation4==0],list1,list2,list3,list4,list5,list6,list7,list8];

In[36]:=
G=(4 \[Pi]^2)/M;
M=m1+m2+m3+m4;
m1=1;
m2=0.000003003;
m3=0.001;
m4=0;
tmax=15;
\[Delta]=0;

In[44]:=
Clear[x1,x2,x3,x4,y1,y2,y3,y4,z1,z2,z3,z4]
Planets=NDSolve[Flatten[{EquationList,WhenEvent[Sqrt[{x4[t],y4[t],z4[t]}.{x4[t],y4[t],z4[t]}]==40.01,{"StopIntegration",tstop=t}]}],Join[r1,r2,r3,r4],{t,0,tmax},
    MaxSteps->100000,AccuracyGoal->10,WorkingPrecision->20]//Quiet;
Planet1[t_]=r1/.Planets;
Planet2[t_]=r2/.Planets;
Planet3[t_]=r3/.Planets;
Planet4[t_]=r4/.Planets



In[50]:= Planet4[0.000340921065169852074016373895628471032722346106772892051`20]
POSTED BY: Brandyn Wessels
Posted 10 years ago

Hi Brandyn,

It would help if you would post code. But . . . NDSolve returns a rule which can be used to define a solution. NDSolveValue returns the solution itself. In both cases the solution is an interpolating function. It may be that the "value" you are seeing is the range limits that are displayed in the representation of the interpolating function. But the function itself accepts values of the independent variable, with the solution "valid" within the displayed range.You can see this if you execute the code below. (I deleted the output before posting because it would be very messy.)

Kind regards, David

eq = {y''[t] == -1, y[0] == 0, y'[0] == 10};

(* NDSolve returns a rule -- use it to define solution variable *)
solution = y /. NDSolve[eq, y, {t, 0, 20}][[1]]

Plot[solution[t], {t, 0, 20}]

(* NDSolveValue returns a solution -- set variable to it *)
solution = NDSolveValue[eq, y, {t, 0, 20}]

Plot[solution[t], {t, 0, 20}]
POSTED BY: David Keith
Reply to this discussion
Community posts can be styled and formatted using the Markdown syntax.
Reply Preview
Attachments
Remove
or Discard

Group Abstract Group Abstract