| Author |
Comment/Response |
Quincy Browne
|
 07/20/10 2:44pm
Hi everyone,
I'm attempting to use Mathematica 7 to run a physical model simulating the motion of a hair in an oscillating medium.
I've already sucessfully implemented the model using SciLab, so I have verified output that I can compare to for debugging.
Using my Scilab model, I've found that my Mathematica notebook will return correct values all the way until NDSolve handles them. So, I've narrowed my problem down to NDSolve. Specifically, the output of NDSolve does not match the output of SciLab's ODE solver. I've tried using Method->"ExplicitRungeKutta" but this fails due to the stiffness of my ODE (something that doesn't come up using an RK solution method in scilab).
Here's my entire code:
T = 30;
Sal = 30;
sal = Sal/1000;
rhow = 999.92293295 + 0.020341179217*T - 0.0061624591598*T^2 +
0.000022614664708*T^3 - 4.6570659168*^-8*T^4;
Drho = 802.00240891*sal - 2.0005183488*(sal*T) +
0.016771024982*sal*T^2 - 0.000030600536746*sal*T^3 -
0.000016132224742*sal^2*T^2;
Rho = Drho + rhow;
meww = 0.000042844324477 +
1/(0.15700386464*(T + 64.99262005)^2 - 91.296496657);
mewA = 1.540913604 + 0.019981117208*T - 0.000095203865864*T^2;
mewB = 7.9739318223 - 0.075614568881*T + 0.00047237011074*T^2;
Mew = meww*(1 + mewA*sal + mewB*sal^2);
L = 0.005;
f = 100;
d = 7.*^-6;
Rhoh = 1100;
R = 0;
S = 4/10^12;
Uo = 0.005;
v = Mew/Rho;
w = 2*Pi*f;
s = (d/4)*((2*Pi*f)/v)^0.5;
g = 0.577 + Log[s];
G = -g/(g^2 + Pi^2/16);
BETA = (w/(2*v))^0.5;
Ih = ((Pi*Rhoh*d^2)/48)*(4*L^3 + (3/4)*d^2*L);
Irho = (Pi*Rho*d^2*L^3)/12;
Imew = -(Pi^2*Mew*G*L^3)/(3*g*w);
Rmew = (4/3)*Pi*Mew*G*L^3;
It = Ih + Irho + Imew;
Rt = Rmew + R;
t = 1;
tmin = 0;
tmax = 0.4;
tint = tmax/10;
tnext = t - tint;
tprev = t + tint;
Td = NIntegrate[
y*Uo*(-Cos[w*t] + Cos[w*t - BETA*y]*Exp[(-BETA)*y]), {y, 0, L}]
Tvm = NIntegrate[
y*Uo*(w*Sin[w*t] - w*Sin[w*t - BETA*y]*Exp[(-BETA)*y]), {y, 0, L}]
c = 4*Pi*Mew*G*Td + ((Pi*Rho*d^2)/4 - (Pi^2*Mew*G)/(g*w))*Tvm;
c = c/It
a = Rt/It
b = S/It
x = NDSolve[{Derivative[2][sol][p] == (-a)*Derivative[1][sol][p] -
b*sol[p] + c, Derivative[1][sol][0] == sol[0] == 0}, {sol}, {p,
t, t}]
solutions = {sol[t], Derivative[1][sol][t], Derivative[2][sol][t]} /.x
Thanks in advance!
URL: , |
|
