MathGroup Archive 2009

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

Search the Archive

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



  • Prev by Date: Re: Re: Re: Select text with keyboard?
  • Next by Date: Re: Logarithmic Axes for ListPlot3D
  • Previous by thread: Mathematica Special Interest Group (Washington DC Area)
  • Next by thread: Re: why the inputs are not ODEs in NDSolve function