Student Support Forum: 'NDSolve Problems' topicStudent Support Forum > General > Archives > "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