MathGroup Archive 2008

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

Search the Archive

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

> ***************************

> 




  • Prev by Date: convex FindMinimum for a large vector
  • Next by Date: Re: C++ const and mathlink
  • Previous by thread: NDSolve[] and Differential Equations: Problem solving two similar
  • Next by thread: multiobjective methods in methamatica