Re: NDSolve[] and Differential Equations: Problem solving two similar
- To: mathgroup at smc.vnet.net
- Subject: [mg85707] Re: NDSolve[] and Differential Equations: Problem solving two similar
- From: dh <dh at metrohm.ch>
- Date: Tue, 19 Feb 2008 07:07:52 -0500 (EST)
- References: <fpdul8$r3l$1@smc.vnet.net>
Hi Gopinath,
in your equations there is a function Phi that is not included in the
list: bdefmitr of functions to solve for.
hope this helps, Daniel
Gopinath Venkatesan wrote:
> 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
> ***************************
>