Mathematica 9 is now available
Services & Resources / Wolfram Forums
-----
 /
MathGroup Archive
2004
*January
*February
*March
*April
*May
*June
*July
*August
*September
*October
*November
*December
*Archive Index
*Ask about this page
*Print this page
*Give us feedback
*Sign up for the Wolfram Insider

MathGroup Archive 2004

[Date Index] [Thread Index] [Author Index]

Search the Archive

ND Solver problem

  • To: mathgroup at smc.vnet.net
  • Subject: [mg46787] ND Solver problem
  • From: "Joseph Lopez" <jo69erbaby at hotmail.com>
  • Date: Mon, 8 Mar 2004 04:10:22 -0500 (EST)
  • Sender: owner-wri-mathgroup at wolfram.com

Hello Steve,

I have been using a program i have written which I am using for a project. 
The program differentiates some equations to obtain second order 
differential equations (ODEs). These equations are then integrated 
numerically using NDsolver to get the results that I want. The program below 
works fine. However when I increase the variables in my initial equations 
and then differentiate them the resultant ODEs are very large. When I 
integrate these using NDsolve an error appears saying that there is a 
singularity in my equations. I was wondering is there a limit to how large 
my ODEs can be without errors occurring from the NDsolver. I am very eager 
to know if there is. My initial program came up with the same error but then 
I simplified the equations and it worked fine. The equations in my other 
program are still quite large after I simplify them however. I would be very 
grateful if you could please give me some feedback on this matter.
The first program below is my initial program the other after it is my 
second program that has more variables.

Thank you

Joseph. Lopez

first Program:

Clear[eqn1,eqn2,eqn3,MEarth,G,e,Rp,R0,\[Theta]start,\[Mu],a,p,VR,R0,orbvel,
    TOrbit,time,ans,R1,M];

\!\(\(cosB = \((R1^2 + R1^2 - R1^2)\)/\((2*R1*R1)\);\)\[IndentingNewLine]
  \(Mt = M + M + M;\)\[IndentingNewLine]
  \(x = \((M*R1 + M*R1*cosB)\)/Mt;\)\[IndentingNewLine]
  \(y = \((M*\@\((R1^2 - \((R1*cosB)\)^2)\))\)/Mt;\)\[IndentingNewLine]
  \(u = \@\((R1^2 - \((R1*cosB)\)^2)\) - y;\)\[IndentingNewLine]
  \(C1 = \@\(\((R1*cosB - x)\)^2 + \((u)\)^2\);\)\[IndentingNewLine]
  \(C2 = \@\(y^2 + x^2\);\)\[IndentingNewLine]
  \(C3 = \@\(\((R1 - x)\)^2 + y^2\);\)\)

\[Phi]2=ArcCos[(C2^2+C1^2-R1^2)/(2*C2*C1)];
\[Phi]3=ArcCos[(C2^2+C3^2-R1^2)/(2*C2*C3)];
x1=R[t]*Cos[\[Theta][t]]+C1*Cos[\[Phi]1[t]+\[Theta][t]];
y1=R[t]*Sin[\[Theta][t]]+C1*Sin[\[Phi]1[t]+\[Theta][t]];
x2=R[t]*Cos[\[Theta][t]]+C2*Cos[\[Phi]2+\[Phi]1[t]+\[Theta][t]];
y2=R[t]*Sin[\[Theta][t]]+C2*Sin[\[Phi]2+\[Phi]1[t]+\[Theta][t]];
x3=R[t]*Cos[\[Theta][t]]+C3*Cos[\[Phi]3+\[Phi]2+\[Phi]1[t]+\[Theta][t]];
y3=R[t]*Sin[\[Theta][t]]+C3*Sin[\[Phi]3+\[Phi]2+\[Phi]1[t]+\[Theta][t]];

\!\(\(\(\[IndentingNewLine]\)\(\(Ro1 = \@\(x1\^2 + y1\^2\);\)\
\[IndentingNewLine]
  \(Ro2 = \@\(x2\^2 + y2\^2\);\)\[IndentingNewLine]
  \(Ro3 = \@\(x3\^2 + y3\^2\);\)\)\)\)

\!\(\(x1vel = \[PartialD]\_t\ x1;\)\[IndentingNewLine]
  \(y1vel = \[PartialD]\_t\ y1;\)\[IndentingNewLine]
  \(x2vel = \[PartialD]\_t\ x2;\)\[IndentingNewLine]
  \(y2vel = \[PartialD]\_t\ y2;\)\[IndentingNewLine]
  \(x3vel = \[PartialD]\_t\ x3;\)\[IndentingNewLine]
  \(y3vel = \[PartialD]\_t\ y3;\)\)

T=1/2*M*(x1vel^2+y1vel^2)+1/2*M*(x2vel^2+y2vel^2)+1/2*M*(x3vel^2+y3vel^2);

\!\(\*
  RowBox[{
    RowBox[{
      RowBox[{"T", "=",
        RowBox[{"(",
          RowBox[{\(1\/2\), " ", "M", " ",
            RowBox[{"(",
              RowBox[{
                RowBox[{"3", " ",
                  SuperscriptBox[
                    RowBox[{
                      SuperscriptBox["R", "\[Prime]",
                        MultilineFunction->None], "[", "t", "]"}], "2"]}],
                "+",
                RowBox[{\((R1\^2 + 3\ R[t]\^2)\), " ",
                  SuperscriptBox[
                    RowBox[{
                      SuperscriptBox["\[Theta]", "\[Prime]",
                        MultilineFunction->None], "[", "t", "]"}], "2"]}],
                "+",
                RowBox[{"2", " ", \(R1\^2\), " ",
                  RowBox[{
                    SuperscriptBox["\[Theta]", "\[Prime]",
                      MultilineFunction->None], "[", "t", "]"}], " ",
                  RowBox[{
                    SuperscriptBox["\[Phi]1", "\[Prime]",
                      MultilineFunction->None], "[", "t", "]"}]}], "+",
                RowBox[{\(R1\^2\), " ",
                  SuperscriptBox[
                    RowBox[{
                      SuperscriptBox["\[Phi]1", "\[Prime]",
                        MultilineFunction->None], "[", "t", "]"}], "2"]}]}],
              ")"}]}], ")"}]}], ";"}], "\[IndentingNewLine]"}]\)

V=-(G*MEarth)*(M/(Ro1)+M/(Ro2)+M/(Ro3));

\!\(\(\(V = \(-G\)\ MEarth\ \((M\ \((1\/\@\(R1\^2\/3 + \(2\ \@R1\^2\ Cos[\
\[Phi]1[t]]\ R[t]\)\/\@3 + R[t]\^2\) +
                1/\((\[Sqrt]\((R1\^2\/3 + R[t]\^2 -
                          1\/3\ \@R1\^2\ R[
                              t]\ \((\@3\ Cos[\[Phi]1[t]] -
                                3\ Sin[\[Phi]1[t]])\))\))\) +
                1/\((\[Sqrt]\((R1\^2\/3 + R[t]\^2 -
                          1\/3\ \@R1\^2\ R[
                              t]\ \((\@3\ Cos[\[Phi]1[t]] +
                                3\ Sin[\[Phi]1[
                                      t]])\))\))\))\))\);\)\(\
\[IndentingNewLine]\)
  \)\)

Lagrange=T-V;

\!\(\(eqn1 = \[PartialD]\_t\ \((\[PartialD]\_\(\(\[Theta]'\)[t]\)Lagrange)\) 
\
- \[PartialD]\_\(\[Theta][t]\)Lagrange;\)\)

\!\(\(eqn2 = \[PartialD]\_t\ \((\[PartialD]\_\(\(R'\)[t]\)Lagrange)\) - \
\[PartialD]\_\(R[t]\)Lagrange;\)\)

\!\(\(eqn3 = \[PartialD]\_t\ \((\[PartialD]\_\(\(\[Phi]1'\)[t]\)Lagrange)\) 
- \
\[PartialD]\_\(\[Phi]1[t]\)Lagrange;\)\)

\!\(\({MEarth, G, e, Rp, \[Theta]start, R1, M} = {5.97*10^24,
        6.672*10^\((\(-11\))\), 0, 7000*10^3, 0, 1*10^3,
        100};\)\[IndentingNewLine]
  \(\[Mu] = G*MEarth;\)\[IndentingNewLine]
  \(a = Rp*\((1 + e)\)/\((1 - e^2)\);\)\[IndentingNewLine]
  \(p = Rp*\((1 + e)\);\)\[IndentingNewLine]
  \(VR = \@\(\[Mu]/p\)*e*Sin[\[Theta]start];\)\[IndentingNewLine]
  \(R0 = p/\((1 + e*Cos[\[Theta]start])\);\)\[IndentingNewLine]
  \(orbvel = \@\(\[Mu]/p\)*\((1 + e*Cos[\[Theta]start])\)/
          R0;\)\[IndentingNewLine]
  \(TOrbit = 2*\[Pi]/\@\(\[Mu]/\((a^3)\)\);\)\[IndentingNewLine]
  \(time = 16*TOrbit;\)\)

ans=NDSolve[{eqn1\[Equal]0,eqn2\[Equal]0,
      eqn3\[Equal]0,\[Theta]'[0]\[Equal]
        orbvel,\[Theta][0]\[Equal]\[Theta]start,R'[0]\[Equal]0,
      R[0]\[Equal]R0,\[Phi]1'[0]\[Equal]0,\[Phi]1[0]\[Equal]\[Pi]/
          4},{\[Theta],R,\[Phi]1},{t,0,time},MaxSteps\[Rule]Infinity]
Plot[Evaluate[\[Theta][t]/.ans],{t,0,time},
  AxesLabel\[Rule]{"Time[sec]","\[Theta][t]"},ImageSize\[Rule]350,
  PlotRange\[Rule]All]
Plot[Evaluate[\[Phi]1'[t]/.ans],{t,0,time},
    AxesLabel\[Rule]{"Time[sec]","\[Phi]'[t]"},ImageSize\[Rule]350,
    PlotRange\[Rule]All]Plot[Evaluate[\[Phi]1[t]/.ans],{t,0,time},
    AxesLabel\[Rule]{"Time[sec]","\[Phi][t]"},ImageSize\[Rule]350,
    PlotRange\[Rule]All]Plot[Evaluate[R[t]/.ans],{t,0,time},
    AxesLabel\[Rule]{"Time[sec]","R[t]"},ImageSize\[Rule]350,
    PlotRange\[Rule]All]

Second Program:

Clear[eqn1,eqn2,eqn3,MEarth,G,e,Rp,\[Theta]start,\[Mu],a,p,VR,R0,orbvel,
    TOrbit,time,ans,R1,R2,R3,M,Mt];

\!\(\(cosB = \((R1^2 + R2^2 - R3^2)\)/\((2*R1*R2)\);\)\n
  \(Mt = M + M + M;\)\n
  \(x = \((M*R2 + M*R1*cosB)\)/Mt;\)\n
  \(y = \((M*\@\((R1^2 - \((R1*cosB)\)^2)\))\)/Mt;\)\n
  \(u = \@\((R1^2 - \((R1*cosB)\)^2)\) - y;\)\n
  \(C1 = \@\(\((R1*cosB - x)\)^2 + \((u)\)^2\);\)\n
  \(C2 = \@\(y^2 + x^2\);\)\n
  \(C3 = \@\(\((R2 - x)\)^2 + y^2\);\)\)

\[Phi]2=ArcCos[(C2^2+C1^2-R1^2)/(2*C2*C1)];
\[Phi]3=ArcCos[(C2^2+C3^2-R2^2)/(2*C2*C3)];
x1=R[t]*Cos[\[Theta][t]]+C1*Cos[\[Phi]1[t]+\[Theta][t]];
y1=R[t]*Sin[\[Theta][t]]+C1*Sin[\[Phi]1[t]+\[Theta][t]];
x2=R[t]*Cos[\[Theta][t]]+C2*Cos[\[Phi]2+\[Phi]1[t]+\[Theta][t]];
y2=R[t]*Sin[\[Theta][t]]+C2*Sin[\[Phi]2+\[Phi]1[t]+\[Theta][t]];
x3=R[t]*Cos[\[Theta][t]]+C3*Cos[\[Phi]3+\[Phi]2+\[Phi]1[t]+\[Theta][t]];
y3=R[t]*Sin[\[Theta][t]]+C3*Sin[\[Phi]3+\[Phi]2+\[Phi]1[t]+\[Theta][t]];

\!\(\(Ro1 = \@\(x1\^2 + y1\^2\);\)\[IndentingNewLine]
  \(Ro2 = \@\(x2\^2 + y2\^2\);\)\[IndentingNewLine]
  \(Ro3 = \@\(x3\^2 + y3\^2\);\)\)

\!\(\(x1vel = \[PartialD]\_t\ x1;\)\[IndentingNewLine]
  \(y1vel = \[PartialD]\_t\ y1;\)\[IndentingNewLine]
  \(x2vel = \[PartialD]\_t\ x2;\)\[IndentingNewLine]
  \(y2vel = \[PartialD]\_t\ y2;\)\[IndentingNewLine]
  \(x3vel = \[PartialD]\_t\ x3;\)\[IndentingNewLine]
  \(y3vel = \[PartialD]\_t\ y3;\)\)

T=1/2*M*(x1vel^2+y1vel^2)+1/2*M*(x2vel^2+y2vel^2)+1/2*M*(x3vel^2+y3vel^2);

V=-(G*MEarth)*(M/(Ro1)+M/(Ro2)+M/(Ro3));

Lagrange=T-V;

\!\(\(eqn1 = \[PartialD]\_t\ \((\[PartialD]\_\(\(\[Theta]'\)[t]\)Lagrange)\) 
\
- \[PartialD]\_\(\[Theta][t]\)Lagrange;\)\)

\!\(\(eqn2 = \[PartialD]\_t\ \((\[PartialD]\_\(\(R'\)[t]\)Lagrange)\) - \
\[PartialD]\_\(R[t]\)Lagrange;\)\)

\!\(\(eqn3 = \[PartialD]\_t\ \((\[PartialD]\_\(\(\[Phi]1'\)[t]\)Lagrange)\) 
- \
\[PartialD]\_\(\[Phi]1[t]\)Lagrange;\)\)

\!\(\({MEarth, G, e, Rp, \[Theta]start, R1, R2, R3, M} = {5.97*10^24,
        6.672*10^\((\(-11\))\), 0.2, 7000*10^3, 0, 100, 100, 100,
        100};\)\[IndentingNewLine]
  \(\[Mu] = G*MEarth;\)\[IndentingNewLine]
  \(a = Rp*\((1 + e)\)/\((1 - e^2)\);\)\[IndentingNewLine]
  \(p = Rp*\((1 + e)\);\)\[IndentingNewLine]
  \(VR = \@\(\[Mu]/p\)*e*Sin[\[Theta]start];\)\[IndentingNewLine]
  \(R0 = p/\((1 + e*Cos[\[Theta]start])\);\)\[IndentingNewLine]
  \(orbvel = \@\(\[Mu]/p\)*\((1 + e*Cos[\[Theta]start])\)/
          R0;\)\[IndentingNewLine]
  \(TOrbit = 2*\[Pi]/\@\(\[Mu]/\((a^3)\)\);\)\[IndentingNewLine]
  \(time = 4*TOrbit;\)\)

ans=NDSolve[{eqn1\[Equal]0,eqn2\[Equal]0,
      eqn3\[Equal]0,\[Theta]'[0]\[Equal]
        orbvel,\[Theta][0]\[Equal]\[Theta]start,R'[0]\[Equal]0,
      R[0]\[Equal]R0,\[Phi]1'[0]\[Equal]0,\[Phi]1[0]\[Equal]0},{\[Theta],
      R,\[Phi]1},{t,0,time},MaxSteps\[Rule]Infinity]
Plot[Evaluate[\[Theta][t]/.ans],{t,0,time},
  AxesLabel\[Rule]{"Time[sec]","\[Phi][t]"},ImageSize\[Rule]350,
  PlotRange\[Rule]All]

_________________________________________________________________
Sign-up for a FREE BT Broadband connection today! 
http://www.msn.co.uk/specials/btbroadband


  • Prev by Date: RE: Shared variables between notebooks?
  • Next by Date: How to install MultipleLogListPlot ?
  • Previous by thread: Re: Differences in Random Numbers
  • Next by thread: How to install MultipleLogListPlot ?