NDSolve[] with error dsfun:
- To: mathgroup at smc.vnet.net
- Subject: [mg90752] NDSolve[] with error dsfun:
- From: Gopinath Venkatesan <gopinathv at ou.edu>
- Date: Wed, 23 Jul 2008 05:56:24 -0400 (EDT)
Hello Mathematica Friends,
When I used NDSolve[] to a set of equations (one equation used an interpolating function), I got an error message saying,
NDSolve::dsfun: {InterpolatingFunction[{{0.,1.}},<<3>>,{Automatic}][t]} cannot be used as a function. >>
Error message details show that somehow function is replaced by string or number. But by replacing with a number still means a meaningful ODE (constant term), and expect it to solve with no problem. But with string, yes there is a problem. But I really don't know why a well defined interpolating function has such a problem. The below given code (for Mathematica version 6) takes only less than 30 seconds to run. So please try it.
Please suggest me ways to solve this problem. Thank you.
Regards,
Gopinath
University of Oklahoma
Below is the code:
(* Starts here *)
mperl = 2303;
ag = 9.81;
bs = 3.757;
ht = 2.1;
moi = (bs ht^3)/12;
emod = 2.87*10^9;
spload = 508329.675;
waload = 127082.41875;
load = spload + waload;
Print["The total load is ", load];
tmass = load/ag;
spmass = spload/ag;
wamass = waload/ag;
Print["The sprung load of vehicle (Chassis) is ", spload,
"and sprung mass is ", spmass];
Print["The unsprung load of vehicle (Wheel/Axle) is ", waload,
"and unsprung mass is ", wamass];
Print["The total load of vehicle is ", load, "and the total mass is ",
tmass];
logdampdecbeam = 0.5;
len = 25;
vel = 30; (* default is 4.778 *)
Tfact = len/vel;
\[Beta] = 1/(2 \[Pi]) logdampdecbeam;
deflmid = (2 load len^3)/(\[Pi]^4 emod moi);
Print["The midspan deflection of SS beam loaded by P = ", load,
" is ", deflmid];
brK = load/deflmid;
Print["Total bridge stiffness is given by ", brK];
\[Omega][j_] := (j^2 \[Pi]^2)/len^2 Sqrt[(emod moi)/mperl];
\[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]];
sK = (0.5 \[Omega][
1])^2 spmass; (* sK is equivalent to parameter C of Fryba's and \
change 0.5 to some value to effect a change in srf, srf \[Congruent] \
0.5 *)
wK = (7.071 \[Omega][
1])^2 wamass; (* wK is equivalent to parameter Subscript[K, 1] \
of Fryba's and change 7.071 to effect a change in wrf, wrf \
\[Congruent] 7.071 *)
fsm = 1/(2 \[Pi]) Sqrt[sK/
spmass]; (* The reason fsm is defined like this instead of fsm = \
0.5 \[Omega][1]/(2 \[Pi]) or fsm = srf \[Omega][1]/(2 \[Pi]) is to \
have flexibility to change sK later, srf \[Congruent] 0.5 *)
fwm = 1/(2 \[Pi]) Sqrt[wK/wamass];
Print["The frequency of sprung (chassis) mass is ", fsm];
Print["The frequency of unsprung (wheel-axle) mass is ", fwm];
dsK = 2 0.5 spmass fsm; (* Here 0.5 stands for the value of \
logarithmic decrement of vehicle spring damping, dwraw \[Congruent] \
0.5 *)
Print["The spring constant of suspension system between chassis and \
wheel-axle system is ", sK];
Print["The viscous damping coefficient \!\(\*SubscriptBox[\"C\", \
\"b\"]\) of suspension system between chassis and wheel-axle system \
is ", dsK];
Print["The spring constant for the tire is ", wK];
\[Kappa] = load/(ag mperl len);
Print["The mass ratio of vehicle to beam (Fryba \[Chi]) is ", \
\[Kappa]];
\[Kappa]1 = waload/spload;
Print["The wheel-axle to chassis mass ratio (Fryba's \
\!\(\*SubscriptBox[\"\[Chi]\", \"0\"]\)) is ", \[Kappa]1];
\[Alpha] = vfreq/\[Omega][1];
Print["The speed parameter is ", \[Alpha]];
srf = (2 \[Pi] fsm)/\[Omega][
1]; (* srf is equivalent to Subscript[\[Gamma], 2] of Fryba's and \
if you want to change srf, change 0.5 in sK definition *)
gammasd = (sK deflmid)/load; (* gammasd is equivalent to \!
\*SubsuperscriptBox[\(\[Gamma]\), \(2\), \('\)] of Fryba's which can \
be alternatively used for defining sprung frequency parameter *)
gammaud = (wK deflmid)/load; (* gammaud is equivalent to \!
\*SubsuperscriptBox[\(\[Gamma]\), \(1\), \('\)] of Fryba's which can \
be used for defining unsprung frequency parameter *)
wrf = (2 \[Pi] fwm)/\[Omega][
1]; (* wrf is equivalent to Subscript[\[Gamma], 1] of Fryba's and \
if you want to change wrf, change 7.071 in wK definition *)
Print["Check factor ", wrf^2,
" should be equal to ", ((1 + \[Kappa]1)/(
2 \[Kappa] \[Kappa]1)) gammaud];
Print["Check factor ", srf^2,
" should be equal to ", ((1 + \[Kappa]1)/(2 \[Kappa])) gammasd];
dwraw = dsK/(2 spmass fsm);
Print["srf equals ", srf];
Print["wrf equals ", wrf];
Print["The logarithmic decrement of vehicle spring damping is ",
dwraw];
Subscript[\[Omega], b] = \[Beta] \[Omega][1];
\[Omega]jd[j_] := Sqrt[(\[Omega][j])^2 - Subscript[\[Omega], b]^2];
step = 0;
nmx = 3;
compareMatrices[A_?MatrixQ, B_?MatrixQ] :=
If[Dimensions[A] == Dimensions[B],
And @@ MapThread[Greater, Flatten /@ {A, B}],
Print["the matrices are not of the same dimension"]];
bdefm = Table[v[im][t], {im, 1, nmx}];
bdefm1t = D[bdefm, t];
bdefm2t = D[bdefm, {t, 2}];
icsp1 = Table[bdefm[[im]] /. t -> 0, {im, 1, nmx}];
icsp2 = Table[1/Tfact bdefm1t[[im]] /. t -> 0, {im, 1, nmx}];
spld[im_, t_] := Sin[im \[Pi] t];
modset = Table[\[Pi]^2 \[Alpha]^2 bdefm2t[[im]] +
2 \[Pi]^3 \[Alpha] \[Beta] bdefm1t[[im]] +
im^4 \[Pi]^4 bdefm[[im]] - 96 spld[im, t], {im, 1, nmx}];
sol1 = Table[
NDSolve[{modset[[im]] == 0, icsp1[[im]] == 0, icsp2[[im]] == 0},
bdefm[[im]], {t, 0, 1}], {im, 1, nmx}];
bdefmsol1 = Table[Chop[bdefm[[im]] /. sol1[[im]]], {im, 1, nmx}];
Plot[bdefmsol1[[1]], {t, 0, 1}]
bdefspld[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[bdefspld[0.5], {t, 0, 1}]
bdefspld[x]
bdefspld[0.5]
bdeftot[x_] := Flatten[Table[bdefspld[x], {t, 0, 1, 0.05}]];
bdeftot[0.5]
time = Table[t, {t, 0, 1, 0.05}];
xlis = Table[xs, {xs, 0, 1, 0.1}];
gridxt = {xlis, time};
bdeftotintp = Table[Chop[bdeftot[x]], {x, 0, 1, 0.1}];
funcSPL[x_, y_] :=
ListInterpolation[bdeftotintp, gridxt, InterpolationOrder -> 3][x,
y];
funcSPL[0.5, 0.5]
bdefspld[0.5] /. t -> 0.5
definstxt = Table[Chop[funcSPL[t, t]], {t, 0, 1, 0.05}];
Print["The deflection of beam at the point of action of load is "];
definstxt // MatrixForm
definstxtitr = definstxt;
funcSPLitr[x_, y_] = funcSPL[x, y];
wdefitr[t_] = v1[t];
sdefitr[t_] =
v3[t]; (* No v2[t] is defined in this code, just to be consistent \
with the sprung mass definition with 2 axle load problem in which \
sprung mass is denoted by v3[t] *)
oldeflBeam = Table[0, {i, 0, 1, 0.1}, {j, 0, 1, 0.05}];
errdeflBeam = Table[1, {i, 0, 1, 0.1}, {j, 0, 1, 0.05}];
mxsize = Dimensions[oldeflBeam][[1]];
mtsize = Dimensions[oldeflBeam][[2]];
(* Start While/Do Loop from here *)
While[compareMatrices[
Drop[Abs[errdeflBeam], {1, mxsize, mxsize - 1}, {1}],
Drop[Abs[0.05 oldeflBeam], {1, mxsize, mxsize - 1}, {1}]],
Print["The value of step is ", step];
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}];
icon2 = Table[1/Tfact bdefmitr1t[[im]] /. t -> 0, {im, 1, nmx}];
wdeficon1 = v1[t] /. t -> 0;
wdeficon2 = D[v1[t], t] /. t -> 0;
sdeficon1 = v3[t] /. t -> 0;
sdeficon2 = D[v3[t], t] /. t -> 0;
reactP[t] = (2 \[Kappa] \[Kappa]1 wrf^2)/(
1 + \[Kappa]1) (Chop[wdefitr[t]] - Chop[funcSPLitr[t, t]]);
If[step != 0,
Print["After solution and before setting reaction is ",
reactP[t] // MatrixForm];
(* reactP[t]=Evaluate[reactP[t]]/.x_?Negative-> 0; *)
reactP[t] = If[Evaluate[reactP[t]] < 0, 0];
Print["Positive/modified reaction force at loadpoint is ",
reactP[t] // MatrixForm]];
v3eq = v3''[t] + (\[Pi]^2 srf^2)/\[Alpha]^2 (v3[t] - v1[t]) + (
dwraw srf)/\[Alpha] (v3'[t] - v1'[t]);
v1eq = v1''[
t] - \[Pi]^2/\[Alpha]^2 ((1 + \[Kappa]1)/(
2 \[Kappa] \[Kappa]1)) (1 -
reactP[t]) + (\[Pi]^2 srf^2)/(\[Alpha]^2 \[Kappa]1) (v1[t] -
v3[t]) + (dwraw srf)/(\[Kappa]1 \[Alpha]) (v1'[t] - v3'[t]);
modsetitr =
Table[\[Pi]^2 \[Alpha]^2 bdefmitr2t[[im]] +
2 \[Pi]^3 \[Alpha] \[Beta] bdefmitr1t[[im]] +
im^4 \[Pi]^4 bdefmitr[[im]] -
Evaluate[96 reactP[t] Sin[im \[Pi] t]], {im, 1, nmx}];
Print["Before block"];
solitr1 =
Table[NDSolve[{v3eq == 0, v1eq == 0, modsetitr[[im]] == 0,
icon1[[im]] == 0, icon2[[im]] == 0, wdeficon1 == wdeficon2 == 0,
sdeficon1 == sdeficon2 == 0}, {bdefmitr[[im]], wdefitr[t],
sdefitr[t]}, {t, 0, 1}], {im, 1, nmx}];
bdefmitrsol1 =
Table[Chop[bdefmitr[[im]] /. solitr1[[im]]], {im, 1, nmx}];
bdefsplditr[x_] := \!\(
\*UnderoverscriptBox[\(\[Sum]\), \(k =
1\), \(nmx\)]\((bdefmitrsol1[\([\)\(k\)\(]\)]\ Sin[
k\ \[Pi]\ x])\)\);
bdeftotitr[x_] := Flatten[Table[bdefsplditr[x], {t, 0, 1, 0.05}]];
Print["First mode contribution is given by the following plot"];
Print[Plot[bdefmitrsol1[[1]], {t, 0, 1}]];
Print["Second mode contribution is given by the following plot"];
Print[Plot[bdefmitrsol1[[2]], {t, 0, 1}]];
Print["Third mode contribution is given by the following plot"];
Print[Plot[bdefmitrsol1[[3]], {t, 0, 1}]];
Print["Combined/effective midspan displacement is given by the \
following plot"];
Print[p1 = Plot[bdefsplditr[0.5], {t, 0, 1}]];
Print["After block"];
bdeftotintpitr = Table[Chop[bdeftotitr[x]], {x, 0, 1, 0.1}];
Print["bdeftotintpitr is ", bdeftotintpitr // MatrixForm];
funcSPLitr[x_, y_] =
ListInterpolation[bdeftotintpitr, gridxt, InterpolationOrder -> 3][
x, y];
definstxtitr = Table[Chop[funcSPLitr[t, t]], {t, 0, 1, 0.05}];
wdefitr[t_] = v1[t] /. solitr1[[1]];
sdefitr[t_] = v3[t] /. solitr1[[1]];
Print[Plot[wdefitr[t], {t, 0, 1}]];
Print[Plot[sdefitr[t], {t, 0, 1}]];
errdeflBeam = bdeftotintpitr - oldeflBeam;
oldeflBeam = bdeftotintpitr;
step = step + 1];
Print["Solution converged and errdeflBeam is not greater than 0.05 \
times oldeflBeam"];
Table[solitr1[[im, All]], {im, 1, nmx}] // MatrixForm
Print["Analysis of Solution"];
(* Ends Here *)