why the inputs are not ODEs in NDSolve function
- To: mathgroup at smc.vnet.net
- Subject: [mg102649] why the inputs are not ODEs in NDSolve function
- From: Haibo Min <yshdfeng at gmail.com>
- Date: Mon, 17 Aug 2009 04:32:08 -0400 (EDT)
Dear group,
I am conducting a simulation by mathematica. The dynamics of the system is a
set of ODEs. However, when expressing them in the NDSolve function, it tells
me that the inputs are not ODEs. I have checked it again and again but still
have no idea why it is so. Please note that in the equations, there is a
Sign function, since anything else seems OK, I doubt whether the problem
arises here. I put the code below, can anyone help me to have a look at it?
Thank you very much!
ThreadEqual = {lhs__} == {rhs__} :> Thread[{lhs} == {rhs}];
J = ( {
{4.35, 0, 0},
{0, 4.337, 0},
{0, 0, 3.664}
} );
S[Px_?(Length@# == 3 &)] := ( {
{0, -Px[[3]], Px[[2]]},
{Px[[3]], 0, -Px[[1]]},
{-Px[[2]], Px[[1]], 0}
} );
Simgmoi[x_?(Length@# == 5 &)] := {1/(1 + e^-x[[1]]), 1/(
1 + e^-x[[2]]), 1/(1 + e^-x[[3]]), 1/(1 + e^-x[[4]]), 1/(
1 + e^-x[[5]])};
W[t_] := {{w11[t], w12[t], w13[t]}, {w21[t], w22[t], w23[t]}, {w31[t],
w32[t], w33[t]}, {w41[t], w42[t], w43[t]}, {w51[t], w52[t],
w53[t]}};
V = ( {
{0.2, 0.3, 0.6, 0.1, 0.6},
{0.5, 0.8, 0.6, 0.7, 0.9},
{0.1, 0.6, 0.5, 0.9, 0.8},
{1.2, 0.3, 1.5, 0.3, 0.6},
{0.5, 0.4, 0.2, 0.5, 0.6},
{0.2, 0.5, 0.5, 0.8, 1.6},
{1.9, 0.3, 0.5, 1.6, 1.9},
{0.9, 0.6, 1.6, 1.2, 1.3},
{1.1, 1.5, 0.2, 0.3, 0.6},
{1.2, 1.5, 0.6, 0.9, 0.8}
} );
e[t_] := {e1[t], e2[t], e3[t]};
r[t_] := {r1[t], r2[t], r3[t]};
T = 5000; c = \[Pi]/T;
k = 35; \[Alpha] = 20; \[Beta] = 40;
Omg[t_] :=
1/50 {-5 c Sin[2 c t], 8 c Sin[4 c t],
4 c Sin[2 c t]}; \[CapitalUpsilon][t_] :=
1/500 {-6 c Sin[8 c t], 18 c Sin[8 c t], 14 c Sin[8 c t]};
x[t_] := {1, Omg[t][[1]], Omg[t][[2]], Omg[t][[3]], Omg'[t][[1]],
Omg'[t][[2]], Omg'[t][[3]], Omg''[t][[1]], Omg''[t][[2]],
Omg''[t][[3]]};
init = {W == ( {
{0.2, 0.3, 0.1},
{0.5, 0.7, 0.3},
{0.5, 0.2, 0.1},
{0.3, 0.2, 0.6},
{0.7, 0.5, 0.1}
} ), e[0] == {0.2, 0.3, 0.4}, r[0] == {0.1, 0.4, 0.5}};
Expr = D[S[J.(Omg[t] + e[t])].(Omg[t] + e[t]),
t] /. {e1'[t] -> r1[t] - \[Alpha] e1[t],
e2'[t] -> r2[t] - \[Alpha] e2[t],
e3'[t] ->
r3[t] - \[Alpha] e3[
t]};(*The replacement is supposed to be: e'[t]=r[t]-\[Alpha] \
e[t], but I didn't figure out how to implement it...*)
de1 = J.r'[t] ==
J.Omg''[t] + \[Alpha] J.(r[t] - \[Alpha] e[t]) +
Expr + \[CapitalUpsilon]'[
t] + (W[t])\[Transpose].Simgmoi[V\[Transpose].x[t]] + (k + 1) r[
t] + \[Beta] Sign[e[t]];
de2 = e'[t] == r[t] - \[Alpha] e[t];
de3 = W'[t] == 10 KroneckerProduct[Simgmoi[V\[Transpose].x[t]], e[t]];
s = NDSolve[
Flatten[Join[{de1, de2, de3}, init] //. ThreadEqual], {w11, w12,
w13, w21, w22, w23, w31, w32, w33, w41, w42, w43, w51, w52, w53,
e1, e2, e3, r1, r2, r3}, {t, 0, 100}, AccuracyGoal -> \[Infinity]]
Haibo