MathGroup Archive 2008

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

Search the Archive

NDSolve[] and Differential Equations: Problem solving two similar

  • To: mathgroup at smc.vnet.net
  • Subject: [mg85671] NDSolve[] and Differential Equations: Problem solving two similar
  • From: Gopinath Venkatesan <gopinathv at ou.edu>
  • Date: Tue, 19 Feb 2008 01:51:22 -0500 (EST)

Hello Friends,

I used NDSolve[] to solve two cases of differential equations, and first set solves with no problem. A similar set did not solve. I defined a variable as a function having two individual arguments for the sake of introducing differentials simply by dash. Like instead of f[x,y], I used f[x][y], so we can define f[x]'[y] for y derivative. I posted below my code, which is lengthy, but I am posting it because it is required to make my question understandable. I might have used lengthy methods (algorithms), etc.

If you look at the below code (code posted at the bottom of the message), you will notice that the two NDSolve[] as given below,

NDSolve[{modset[[im]] == 0, ictp1[[im]] == 0, ictp2[[im]] == 0}, 
 bdefm[[im]], {t, 0, 1 + sep/len}]

solves with no problem but a similar NDSolve[] as shown below does not.

NDSolve[{eqn1[[im]] == 0, icon1[[im]] == 0, icon2[[im]] == 0}, 
 bdefmitr[[im]], {t, 0, 1 + sep/len}]

Of course, the second equation deals with a set of few dependent differential equations while the first one deals with only one differential equation.

The initial conditions (2 in the first and 10 in the second case) with respective set of differential equations should solve with no problem.

If any of could help me figure out the problem, I will be greatly relieved. Please give me a hint on solving such connected differential equations or if possible let me know where I made the mistake or sources of examples solving dependent multiple differential equations.

Thank you.

Gopinath Venkatesan
University of Oklahoma

******************************
Mathematica code posted below
******************************

mperl = 2303;
ag = 9.81;
bs = 3.757;
ht = 2.1;
moi = (bs ht^3)/12;
emod = 2.87*10^9;
sep = 8;
pmI = 11515 sep^2; 
spload = 508329.675;
fwload = 28240.5375;
rwload = 28240.5375;
d1 = 0.33 sep;
d2 = sep - d1;
load = spload + fwload + rwload;
Print["The total load is ", load];
tmass = load/ag;
logdampdecbeam = 0.5;

len = 25;
vel = 30; (* default is 4.778 *)
Tfact = len/vel;
mfact = 1/(2 - sep/len);
\[Beta] = 1/(2 \[Pi]) logdampdecbeam;
dwrf1 = 0.5;
dwrf2 = 0.5;

deflmid = (2 load len^3)/(\[Pi]^4 emod moi);

\[Omega][j_] := (j^2 \[Pi]^2)/len^2 Sqrt[(emod moi)/mperl];
q1 = (fwload + spload d2/sep)/load;
q2 = (rwload + spload d1/sep)/load;

deflmid1 = q1 deflmid;
deflmid2 = q2 deflmid;

sK1 = (0.04 \[Omega][1])^2 tmass; 
sK2 = (0.04 \[Omega][1])^2 tmass;
\[CapitalOmega][k_] := (k \[Pi])/Tfact;
vfreq = (\[Pi] vel)/len;
w = (\[Pi] vel)/len; (* Remove this parameter later *)
Print["The natural frequency of the beam is ", \[Omega][1]];
fs1 = 1/(2 \[Pi]) Sqrt[sK1/
   tmass]; (* The reason fs1 is defined like this instead of fs1 = \
0.04 \[Omega][1]/(2 \[Pi]) or fs1 = srf1 \[Omega][1]/(2 \[Pi]) is to \
have flexibility to change sK1 later *)
fs2 = 1/(2 \[Pi]) Sqrt[sK2/tmass];

dsK1 = 2 dwrf1 tmass fs1;
dsK2 = 2 dwrf2 tmass fs2;
wK1 = (0.2 \[Omega][
      1])^2 tmass; (* wK1 is equivalent to parameter Subscript[K, 1] \
of Fryba's and change 0.2 to effect a change in wrf1 *)
wK2 = (0.2 \[Omega][1])^2 tmass;

fw1 = 1/(2 \[Pi]) Sqrt[wK1/tmass];

fw2 = 1/(2 \[Pi]) Sqrt[wK2/tmass];

\[Kappa] = load/(ag mperl len);
\[Kappa]1 = fwload/load;
\[Kappa]2 = rwload/load;
rI = (pmI ag)/(load sep^2);

\[Alpha] = vfreq/\[Omega][1];
Print["The speed parameter is ", \[Alpha]];
srf1 = (2 \[Pi] fs1)/\[Omega][1];
srf2 = (2 \[Pi] fs2)/\[Omega][1];
wrf1 = (2 \[Pi] fw1)/\[Omega][1];
wrf2 = (2 \[Pi] fw2)/\[Omega][1];

Print["srf1 equals ", srf1];
Print["srf2 equals ", srf2];
Print["wrf1 equals ", wrf1];
Print["wrf2 equals ", wrf2];

Subscript[\[Omega], b] = \[Beta] \[Omega][1];
\[Omega]jd[j_] := Sqrt[(\[Omega][j])^2 - Subscript[\[Omega], b]^2];

dwrf1 = dsK1/(2 tmass fs1);
dwrf2 = dsK2/(2 tmass fs2);

step = 0;
nmx = 3;

bdefm = Table[v[im][t], {im, 1, nmx}];

bdefm1t = D[bdefm, t];
bdefm2t = D[bdefm, {t, 2}];

ictp1 = Table[bdefm[[im]] /. t -> 0, {im, 1, nmx}];
ictp2 = Table[1/Tfact bdefm1t[[im]] /. t -> 0, {im, 1, nmx}];


tpld[im_, t_] := 
  If[0 <= t <= sep/len, 2 q1 Sin[im \[Pi] t], 
   If[sep/len < t < 1, 
    2 q1 Sin[im \[Pi] t] + 2 q2 Sin[im \[Pi] (t - sep/len)], 
    If[1 <= t <= 1 + sep/len, 2 q2 Sin[im \[Pi] (t - sep/len)], 0]]];

modset = Table[\[Pi]^2 \[Alpha]^2 bdefm2t[[im]] + 
    2 \[Pi]^3 \[Alpha] \[Beta] bdefm1t[[im]] + 
    im^4 \[Pi]^4 bdefm[[im]] - 48 tpld[im, t], {im, 1, nmx}];

sol1 = Table[
   NDSolve[{modset[[im]] == 0, ictp1[[im]] == 0, ictp2[[im]] == 0}, 
    bdefm[[im]], {t, 0, 1 + sep/len}], {im, 1, nmx}];

bdefmsol1 = Table[Chop[bdefm[[im]] /. sol1[[im]]], {im, 1, nmx}];

Plot[bdefmsol1[[1]], {t, 0, 1 + sep/len}]

bdeftpld[x_] := \!\(
\*UnderoverscriptBox[\(\[Sum]\), \(im = 
     1\), \(nmx\)]\((bdefmsol1[\([\)\(im\)\(]\)]\ Sin[
      im\ \[Pi]\ x])\)\); (* this is the beam displacement solution \
we use to substitute later for iterating the vehicle beam interaction \
*)

p1 = Plot[bdeftpld[0.5], {t, 0, 1 + sep/len}]
bdeftpld[x]
bdeftpld[0.5]

bdeftot[x_] := 
  Flatten[Table[bdeftpld[x], {t, 0, 1 + sep/len, 0.05}]];
bdeftot[0.5]

time = Table[t, {t, 0, 1 + sep/len, 0.05}];
xlis = Table[xs, {xs, 0, 1, 0.1}];
gridxt = {xlis, time};

bdeftotintp = Table[Chop[bdeftot[x]], {x, 0, 1, 0.1}];

func[x_, y_] := 
  ListInterpolation[bdeftotintp, gridxt, InterpolationOrder -> 3][x, 
   y];

func[0.5, 0.5]
bdeftpld[0.5] /. t -> 0.5

definstxt1 = 
  Table[If[0 <= t <= 1, Chop[func[t, t]], 0], {t, 0, 1 + sep/len, 
    0.05}];
definstxt2 = 
  Table[If[sep/len <= t <= 1 + sep/len, Chop[func[t - sep/len, t]], 
    0], {t, 0, 1 + sep/len, 0.05}];
Print["The deflection of beam at locations of contact of load 1 and 2 \
are "];
definstxt1 // MatrixForm
definstxt2 // MatrixForm
definstxt1itr = definstxt1;
definstxt2itr = definstxt2;

wdef1itr[t_] := v1[t];
wdef2itr[t_] := v2[t];

(* Start While/Do Loop from here *)

bdefmitr = Table[vi[im][t], {im, 1, nmx}];
bdefmitr1t = D[bdefmitr, t];
bdefmitr2t = D[bdefmitr, {t, 2}];
Print["bdefmitr is ", bdefmitr];
Print["bdefmitr1t is ", bdefmitr1t];
Print["bdefmitr2t is ", bdefmitr2t];

icon1 = Table[bdefmitr[[im]] /. t -> 0, {im, 1, nmx}];
Print["icon1 is ", icon1 // MatrixForm];
icon2 = Table[1/Tfact bdefmitr1t[[im]] /. t -> 0, {im, 1, nmx}];
Print["icon2 is ", icon2 // MatrixForm];
wdef1icon1 = v1[t] /. t -> 0;
wdef1icon2 = D[v1[t], t] /. t -> 0;
wdef2icon1 = v2[t] /. t -> 0;
wdef2icon2 = D[v2[t], t] /. t -> 0;
sdeficon1 = v3[t] /. t -> 0;
sdeficon2 = D[v3[t], t] /. t -> 0;
phicon1 = \[Phi][t] /. t -> 0;
phicon2 = D[\[Phi][t], t] /. t -> 0;

wnet1[t_] := v3[t] + d1/sep \[Phi][t] - v1[t];
wnet2[t_] := v3[t] - d2/sep \[Phi][t] - v2[t];

Print["Defining the reaction forces"];
reactP1[t_] := (\[Pi]^4 \[Kappa] wrf1^2)/
   24 (Chop[wdef1itr[t]] - Chop[func[t, t]]);

(* See if we are using definstxt that is updating within loop *)
reactP2[t_] := (\[Pi]^4 \[Kappa] wrf2^2)/
   24 (Chop[wdef2itr[t]] - Chop[func[t - sep/len, t]]);
If[step != 0,
  Print["After solution and before setting reaction 1 is ", 
   reactP1[t] // MatrixForm];
  Print["After solution and before setting, reaction 2 is ", 
   reactP2[t] // MatrixForm];
  
  reactP1[t] = reactP1[t] /. x_?Negative -> 0;
  reactP2[t] = reactP2[t] /. x_?Negative -> 0;
  
  Print["reaction force at load1 is ", reactP1 // MatrixForm];
  Print["reaction force at load 2 position is ", 
   reactP2 // MatrixForm]];

tplditr[im_, t_] := 
  If[0 <= t < sep/len, reactP1[t] Sin[im \[Pi] t], 
   If[sep/len <= t <= 1, 
    reactP1[t] Sin[im \[Pi] t] + 
     reactP2[t] Sin[im \[Pi] (t - sep/len)], 
    If[1 < t <= 1 + sep/len, 
     reactP2[t] Sin[im \[Pi] (t - sep/len)]]]];


phieq = -\[Phi]''[t] - (d1 sK1 Tfact^2)/(sep rI tmass) wnet1[t] - (
    d1 dsK1 Tfact)/(sep rI tmass) wnet1'[t] + (d2 sK2 Tfact^2)/(
    sep rI tmass) wnet2[t] + (d2 dsK2 Tfact)/(sep rI tmass)
     wnet2'[t];

v3eq = (1 - \[Kappa]1 - \[Kappa]2)/Tfact^2 \[Alpha]^2 v3''[t] + 
   srf1^2 \[Pi]^2 wnet1'[t] + dwrf1 \[Alpha] srf1 wnet1'[t] + 
   srf2^2 \[Pi]^2 wnet2[t] + dwrf2 \[Alpha] srf2 wnet2'[t];

v1eq = v1''[t] - (
    24 Tfact^2)/(\[Pi]^2 \[Kappa]1 \[Kappa] \[Alpha]^2) (2 q1 + \
(\[Pi]^4 \[Kappa] srf1^2 wnet1[t])/
      24 + (\[Pi]^2 dwrf1 srf1 \[Kappa] \[Alpha] wnet1'[t])/24 - 
      reactP1[t]);

v2eq = v2''[t] - (
    24 Tfact^2)/(\[Pi]^2 \[Kappa]2 \[Kappa] \[Alpha]^2) (2 q2 + \
(\[Pi]^4 \[Kappa] srf2^2 wnet2[t])/
      24 + (\[Pi]^2 dwrf2 srf2 \[Kappa] \[Alpha] wnet2'[t])/24 - 
      reactP2[t]);


modsetitr = 
  Table[\[Pi]^2 \[Alpha]^2 bdefmitr2t[[im]] + 
    2 \[Pi]^3 \[Alpha] \[Beta] bdefmitr1t[[im]] + 
    im^4 \[Pi]^4 bdefmitr[[im]] - 48 tplditr[im, t], {im, 1, nmx}];
modsetitr // MatrixForm

eqn1 = Table[
   Flatten[{wdef1icon1, wdef1icon2, wdef2icon1, wdef2icon2, sdeficon1,
      sdeficon2, phicon1, phicon2, phieq, v3eq, v1eq, v2eq, 
     modsetitr[[im]]}], {im, 1, nmx}];

Print["Solving the set of equations"];
solnitr1 = 
  Table[NDSolve[{eqn1[[im]] == 0, icon1[[im]] == 0, icon2[[im]] == 0},
     bdefmitr[[im]], {t, 0, 1 + sep/len}], {im, 1, nmx}];

sol1itr1[[1]] // MatrixForm

***************************
END HERE
***************************


  • Prev by Date: Re: Re: "Assuming"
  • Next by Date: Re: [functional approach should give] an even faster way to normalize
  • Previous by thread: Re: Multiple data sets with ListPlot and different PointSizes - Mesh
  • Next by thread: Re: NDSolve[] and Differential Equations: Problem solving two similar