About solution of complex system of DEs - to administrator -this is
- To: mathgroup at smc.vnet.net
- Subject: [mg95232] About solution of complex system of DEs - to administrator -this is
- From: rosenm at mail.bg
- Date: Sun, 11 Jan 2009 06:39:31 -0500 (EST)
Dear colleagues, i want to solve the system of 5 DEs, listed below. I
am using Mathematica 7 and It is solved successfully only if
SolveDealyed is set to True .If it is set to False, i don't receive
any solution. I would like to use different integration methods,
especially Runge-Kutta 4 which is not active when SolveDealyed is set
to True. What can i do in this case?
Apply[Clear, Names["Global`*"]];
A01 = {{Cos[q1[t]], -Sin[q1[t]], 0, 0}, {Sin[q1[t]], Cos[q1[t]], 0,
0}, {0, 0, 1, 0}, {0, 0, 0, 1}};
A12 = {{Cos[q2[t]], 0, Sin[q2[t]], 0}, {0, 1, 0, 0}, {-Sin[q2[t]], 0,
Cos[q2[t]], L1}, {0, 0, 0, 1}};
A23 = {{Cos[q3[t]], 0, Sin[q3[t]], L2}, {0, 1, 0, 0}, {-Sin[q3[t]], 0,
Cos[q3[t]], 0}, {0, 0, 0, 1}};
A34 = {{Cos[q4[t]], 0, Sin[q4[t]], L3}, {0, 1, 0, 0}, {-Sin[q4[t]], 0,
Cos[q4[t]], 0}, {0, 0, 0, 1}};
A45 = {{Cos[q5[t]], -Sin[q5[t]], 0, 0}, {Sin[q5[t]], Cos[q5[t]], 0,
0}, {0, 0, 1, 0}, {0, 0, 0, 1}};
A02 = A01.A12;
A03 = A01.A12.A23;
A05 = A01.A12.A23.A34.A45;
r22 = {L2/2, 0, L2/7, 1}; rB2 = {L2, 0, 0, 1};
r33 = {L3/2, 0, L3/7, 1}; rC3 = {L3, 0, 0, 1};
r55 = {L4, 0, 0, 1};
r20 = A02.r22; rB0 = A02.rB2;
r30 = A03.r33; rC0 = A03.rC3;
r50 = A05.r55;
V20X = D[r20[[1]], t];
V20Y = D[r20[[2]], t];
V20Z = D[r20[[3]], t];
V20 = Sqrt[V20X^2 + V20Y^2 + V20Z^2];
V30X = D[r30[[1]], t];
V30Y = D[r30[[2]], t];
V30Z = D[r30[[3]], t];
V30 = Sqrt[V30X^2 + V30Y^2 + V30Z^2];
V50X = D[r50[[1]], t];
V50Y = D[r50[[2]], t];
V50Z = D[r50[[3]], t];
V50 = Sqrt[V50X^2 + V50Y^2 + V50Z^2];
Omega011 = {0, 0, q1'[t], 0};
Omega122 = {0, q2'[t], 0, 0};
Omega233 = {0, q3'[t], 0, 0};
Omega010 = A01.Omega011;
Omega120 = A02.Omega122;
Omega230 = A03.Omega233;
OmegaA010 = Omega010;
OmegaA120 = Omega010 + Omega120;
OmegaA230 = Omega120 + Omega230;
OmegaA011 = Inverse[A01].OmegaA010;
OmegaA122 = Inverse[A02].OmegaA120;
OmegaA233 = Inverse[A03].OmegaA230;
\!\(\*
TagBox[
RowBox[{
TagBox[
RowBox[{"J1", "=",
RowBox[{"(", "", GridBox[{
{"0", "0", "0", "0"},
{"0", "0", "0", "0"},
{"0", "0", "J1zz", "0"},
{"0", "0", "0", "0"}
}], "", ")"}]}],
Function[BoxForm`e$,
MatrixForm[BoxForm`e$]]], ";",
RowBox[{"J2", "=",
RowBox[{"(", "", GridBox[{
{"J2xx", "J2yx", "J2zx", "0"},
{"J2xy", "J2yy", "J2zy", "0"},
{"J2xz", "J2yz", "J2zz", "0"},
{"0", "0", "0", "0"}
}], "", ")"}]}]}],
Function[BoxForm`e$,
MatrixForm[BoxForm`e$]]]\); \!\(\*
TagBox[
RowBox[{"J3", "=",
RowBox[{"(", "", GridBox[{
{"J3xx", "J3yx", "J3zx", "0"},
{"J3xy", "J3yy", "J3zy", "0"},
{"J3xz", "J3yz", "J3zz", "0"},
{"0", "0", "0", "0"}
}], "", ")"}]}],
Function[BoxForm`e$,
MatrixForm[BoxForm`e$]]]\);
J1zz = 0.92;
J2xx = 2.04; J2yy = 34.8; J2zz = 34.4; J2xy = J2yx = 0; J2xz =
J2zx = 0.13; J2yz = J2zy = 0;
J3xx = 0.43; J3yy = 39.6; J3zz = 39.6; J3xy = J3yx = 0; J3xz =
J3zx = 0.15; J3yz = J3zy = 0;
Mtheta2 =
Fhc*H1*H2*
Sin[q2[t] + 1.46]/Sqrt[H1^2 + H2^2 - 2*H1*H2*Cos[q2[t] + 1.46]];
KR = 1/2*OmegaA011.J1.OmegaA011 + 1/2*OmegaA122.J2.OmegaA122 +
1/2*OmegaA233.J3.OmegaA233;
KT = 1/2*m2*V20^2 + 1/2*m3*V30^2 + 1/2*m5*V50^2;
K = KT + KR;
P = m2*g*r20[[3]] + m3*g*r30[[3]] + m5*g*r50[[3]] + c1*(q1[t]^2)/2 +
c2*(-q2[t] + q2n)^2/2 + c3*(-q3[t] + q3n)^2/2;
F = K - P;
eq1 = D[D[F, q1'[t]], t] - D[F, q1[t]] + b1*q1'[t];
eq2 = D[D[F, q2'[t]], t] - D[F, q2[t]] + b2*q2'[t];
eq3 = D[D[F, q3'[t]], t] - D[F, q3[t]] + b3*q3'[t];
eq4 = D[D[F, q4'[t]], t] - D[F, q4[t]];
eq5 = D[D[F, q5'[t]], t] - D[F, q5[t]];
H1 = 1.725;
H2 = 0.53;
Fhc = 42500;
Te = 10;
m2 = 120.; m3 = 35.; m5 = 300; L1 = 0.6; L2 = 1.7; L3 = 1.65; L4 = 2;
g = 9.81; c1 = 1.2*10^6; c2 = 1.2*10^6; c3 = 1.2*10^6; b1 =
1.4*10^5; b2 = 1.4*10^5; b3 = 1.4*10^5;
q2n = -40*Pi/180; q3n = 60*Pi/180;
ShowStatus[status_] :=
LinkWrite[$ParentLink,
SetNotebookStatusLine[FrontEnd`EvaluationNotebook[],
ToString[status]]];
Timing[sol =
NDSolve[{eq1 == 0, eq2 == 0, eq3 == 0, eq4 == 0, eq5 == 0,
q1[0] == 0., q2[0] == q2n, q3[0] == q3n,
q4[0] == Pi/3 - q2n - q3n, q5[0] == 1, q1'[0] == 0.,
q2'[0] == 0., q3'[0] == 0, q4'[0] == 0, q5'[0] == 0}, {q1, q2,
q3, q4, q5}, {t, 0, Te}, SolveDelayed -> True,
StepMonitor :> ShowStatus[t]]];
{Plot[Evaluate[q1[t]] /. sol, {t, 0, Te}, AxesLabel -> {t, q1}],
Plot[Evaluate[q1'[t]] /. sol, {t, 0, Te}, AxesLabel -> {t, q1'}],
Plot[Evaluate[q1''[t]] /. sol, {t, 0, Te}, AxesLabel -> {t, q1''}]}