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