Mathematica 9 is now available
Student Support Forum
-----
Student Support Forum: 'NDSolve Problems' topicStudent Support Forum > General > "NDSolve Problems"

Next Comment >Help | Reply To Topic
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: ,

Subject (listing for 'NDSolve Problems')
Author Date Posted
NDSolve Problems Quincy Browne 07/20/10 2:44pm
Re: NDSolve Problems Quincy Browne 07/29/10 6:45pm
Next Comment >Help | Reply To Topic