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 *)