Re: Re: NDSolve[] with nested If[] and Piecewise[] usage
- To: mathgroup at smc.vnet.net
- Subject: [mg91501] Re: [mg91481] Re: NDSolve[] with nested If[] and Piecewise[] usage
- From: DrMajorBob <drmajorbob at att.net>
- Date: Sun, 24 Aug 2008 07:06:08 -0400 (EDT)
- References: <12213370.1219474340300.JavaMail.root@m08>
- Reply-to: drmajorbob at longhorns.com
The first NDSolve has no errors.
The second NDSolve output clearly reveals that tplditr, v1, v2, and v3 are
not defined. (Remove the semicolon to see it.) You're solving ONLY for
vi[1][t], so everything else must be defined and numeric.
The third Table of NDSolves (without the semicolon) outputs mentions of
tplditrPiece, v1, v2, and v3 -- still undefined -- and you're solving only
for vi[i][t].
Et cetera.
DSolve has undefined tplditr, v1, v2, v3, et cetera to deal with, and
you're solving only for vi[1][t].
The last statement again has NDSolve choking on undefined functions.
Send a notebook next time. I just spent half an hour eliminating garbage
caused by copy/paste.
And, if you want anyone to read and understand your code, get rid of 95%
of those Evaluates and anything else that's not needed, too.
Bobby
On Sat, 23 Aug 2008 00:42:25 -0500, Gopinath Venkatesan <gopinathv at ou.edu>
wrote:
> 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 *)
>
>
--
DrMajorBob at longhorns.com