Mathematica 9 is now available
Services & Resources / Wolfram Forums / MathGroup Archive
-----

MathGroup Archive 2008

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

Search the Archive

Re: NDSolve[] with nested If[] and Piecewise[] usage

  • To: mathgroup at smc.vnet.net
  • Subject: [mg91481] Re: NDSolve[] with nested If[] and Piecewise[] usage
  • From: Gopinath Venkatesan <gopinathv at ou.edu>
  • Date: Sat, 23 Aug 2008 01:42:25 -0400 (EDT)
  • Organization: The Math Forum

Jean, DrMajorBob

Jean: When you pointed out some of the variables are not defined - I was over confident that there was no mistakes. I  apologize for ignoring your comments.

Yes, after DrMajorBob's reply, I checked back and corrected  those mistakes by including wdef1loop and wdef2loop. Actually we don't need wdef1loop or wdef1itr (and so on), since wdef1itr[t_]:=v1[t] means wdef1itr and v1 are same but I did it this way because I was working on another model (differential quadrature) and the code is almost the same so modified where-ever needed and so those extra definitions. So wdef1itr is same as wdef1loop which is same as v1, and likewise other repetitions are also there.

Thanks for helping me. I hope you will offer more suggestions after reading this post. Thank you.

Gopinath
University of Oklahoma

Here is the corrected version:

(* Code starts here *)

mperl = 2303;
ag = 9.81;
bs = 3.757;
ht = 2.1;
moi = (bs ht^3)/12;
emod = 2.87*10^9;
sep = 8; (* sep is the axle spacing or parameter D of Fryba's *)
pmI = 11515 sep^2; (* pmI is equivalent to parameter I (polar moment \
of inertia) of Fryba's *)
spload = 508329.675;
fwload = 28240.5375;
rwload = 28240.5375;
d1 = 0.33 sep;
d2 = sep - d1;
load = spload + fwload + rwload;
Print["The total load of vehicle 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; (* replace logdampdecbeam with single parameter \
\[Beta] later *)
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; (* sK1 is equivalent to parameter Subscript[C, 1] \
of Fryba's and change 0.04 to some value to effect a change in srf1 *)


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];

Print["The frequency of sprung mass w.r.t to front wheel is ", fs1];
Print["The frequency of sprung mass w.r.t to rear wheel is ", fs2];

Print["The frequency of unsprung front wheel mass is ", fw1];
Print["The frequency of unsprung rear wheel mass is ", fw2];


\[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]; (* srf1 is equivalent to Subscript[\[Beta], 1] of Fryba's and \
if you want to change srf1, change 0.04 in sK1 definition *)
srf2 = (2 \[Pi] fs2)/\[Omega][1];
wrf1 = (2 \[Pi] fw1)/\[Omega][
   1]; (* wrf1 is equivalent to Subscript[\[Gamma], 1] of Fryba's and \
if you want to change wrf1, change 0.2 in wK1 definition *)
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[Evaluate[0 <= t <= sep/len], Evaluate[2 q1 Sin[im \[Pi] t]], 
   Evaluate[
    If[Evaluate[sep/len < t < 1], 
     Evaluate[
      2 q1 Sin[im \[Pi] t] + 2 q2 Sin[im \[Pi] (t - sep/len)]], 
     Evaluate[
      If[Evaluate[1 <= t <= 1 + sep/len], 
       Evaluate[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}];
bdefm1tsol1 = 
  Table[Chop[D[bdefmsol1[[im]], t] /. sol1[[im]]], {im, 1, nmx}];

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}]

bdeftot[x_] := 
  Table[-bdeftpld[x], {t, 0, 1 + sep/len, (len + sep)/(20 len)}];

time = Table[t, {t, 0, 1 + sep/len, (len + sep)/(20 len)}];
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];

pbtpldint = 
 Plot[-func[1/2, t], {t, 0, 1 + sep/len}, PlotRange -> All, 
  PlotStyle -> {Red}]

definstxt1 = 
  Table[If[0 <= t <= 1, Chop[func[t, t]], 0], {t, 0, 1 + sep/len, (
    len + sep)/(20 len)}];
definstxt2 = 
  Table[If[sep/len <= t <= 1 + sep/len, Chop[func[t - sep/len, t]], 
    0], {t, 0, 1 + sep/len, (len + sep)/(20 len)}];
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];
wdef1loop[t_] := wdef1itr[t];
wdef2loop[t_] := wdef2itr[t];

(* Start While/Do Loop from here *)


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

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[wdef1loop[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[wdef2loop[t]] - Chop[func[t - sep/len, t]]);

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

tplditrPiece[im_?NumericQ, t_?NumericQ] := 
  Piecewise[{{reactP1[t] Sin[im \[Pi] t], 
     0 <= t < sep/len}, {reactP1[t] Sin[im \[Pi] t] + 
      reactP2[t] Sin[im \[Pi] (t - sep/len)], 
     sep/len <= t <= 1}, {reactP2[t] Sin[im \[Pi] (t - sep/len)], 
     1 < t <= 1 + 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]);

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

modsetitrPiece = 
  Table[\[Pi]^2 \[Alpha]^2 bdefmitr2t[[im]] + 
    2 \[Pi]^3 \[Alpha] \[Beta] bdefmitr1t[[im]] + 
    im^4 \[Pi]^4 bdefmitr[[im]] - 
    Evaluate[48 tplditrPiece[im, t]], {im, 1, nmx}];

Print["The modsetitr equation is ", modsetitrIfloop // MatrixForm];

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

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

Print["eqn1 is given by ", eqn1 // MatrixForm];
Print["eqn1[[1]] is given by ", eqn1[[1]] // MatrixForm];
Print["eqn1 is given by ", eqn2 // MatrixForm];

Print["Solving the set of equations"];
solnitr1 = 
  Table[NDSolve[{Flatten[eqn1[[im]] == 0], icon1[[im]] == 0, 
     icon2[[im]] == 0}, bdefmitr[[im]], {t, 0, 1 + sep/len}], {im, 1, 
    nmx}];
solnitr2 = 
  Table[NDSolve[{Flatten[eqn2[[im]] == 0], icon1[[im]] == 0, 
     icon2[[im]] == 0}, bdefmitr[[im]], {t, 0, 1 + sep/len}], {im, 1, 
    nmx}];
Print["After Solving the set of equation"];
Print["The solution to first set of Equations - eqn1 is given by ", 
  solnitr1 // MatrixForm];
Print["The solution to second set of Equations - eqn2 is given by ", 
  solnitr2 // MatrixForm];

Print["Other tries"];
eqn3 = Flatten[{icon1[[1]] == 0, icon2[[1]] == 0, wdef1icon1 == 0, 
    wdef1icon2 == 0, wdef2icon1 == 0, wdef2icon2 == 0, sdeficon1 == 0,
     sdeficon2 == 0, phicon1 == 0, phicon2 == 0, phieq == 0, 
    v3eq == 0, v1eq == 0, v2eq == 0, 
    Evaluate[modsetitrIfloop[[1]]] == 0}];
solnitr3 = DSolve[eqn3, {bdefmitr[[1]]}, t]
NDSolve[{eqn1[[1]] == 0, icon1[[1]] == 0, icon2[[1]] == 0}, 
 bdefmitr[[1]], {t, 0, 1}]

(* Code ends here *)


  • Prev by Date: Re: Non Default stylesheet format and syntax coloring
  • Next by Date: four dimensional plot
  • Previous by thread: Re: Partial differential equation with evolving boundary conditions
  • Next by thread: Re: Re: NDSolve[] with nested If[] and Piecewise[] usage