Message Boards Message Boards

Need help with NDSolve

Posted 10 years ago

I need to solve two state variables using ODEs. I have set up the problem using NDSolve and proper boundary conditions. I am getting errors while Mathematica is trying to solve the set of 2 ODEs. I would appreciate any help in figuring out the problem.

R = 518.3; k = 1.304; tp = 10.; z =.;
V1 = .34; p1 = 25 10^6; T1 = 300;
\[Mu] = .002; Lpipe = 10.; dpipe = 3*25.4 10^-3; Rf = (
 128 \[Mu] Lpipe)/(\[Pi] dpipe^4);
V = V1; \[Rho]a = 1.225; Kf = 1.8 10^9;
c = Sqrt[Kf/\[Rho]f]; \[Rho]f = 872; Ca = V/(\[Rho]a c^2);
Lt = 3*25.4 10^-3; dt = 
 3*25.4 10^-3; At = \[Pi]/4 dt^4; Ithroat = (\[Rho]f Lt)/At;

f = 5; n =.; t =.; \[Chi] = 1/f; \[Omega] = 2 \[Pi] f; r = 1.; L = 
 5 r; \[Theta] = \[Omega] t;

\[Theta]2 = \[Theta] + 2 \[Pi]/3.;
\[Theta]3 = \[Theta]2 + 2 \[Pi]/3;

Clear[y1, y2, s1, P1, P2, P3];
P1[\[Theta]_] := -r Sin[\[Theta]] - (r^2 Sin[\[Theta]] Cos[\[Theta]])/
    Sqrt[L^2 - r^2 (Sin[\[Theta]])^2];
P2[\[Theta]_] := P1[\[Theta] + 2 Pi/3];
P3[\[Theta]_] := P1[\[Theta] + 4 Pi/3];

fsum[\[Theta]_] := 
  Max[P1[\[Theta]], 0] + Max[P2[\[Theta]], 0] + Max[P3[\[Theta]], 0];

ClearAll[Qin];

Qin[t_?NumericQ] := fsum[\[Theta]];


(*----- system eqns------*)

(* ------ Eqn1: y1 = pf ------*)

Eqn1 = y1'[t] + y1[t] Rf/Ithroat + y2[t] 1 /Ca -  Rf Qin[t];

(* ------ Eqn 2: y2 = qf ------*)

Eqn2 = y2'[t] - y1[t]  (1/Ithroat);


s1 = NDSolve[{Eqn1 ==  0, Eqn2 == 0, y1[0] == 10 10^6, 
    y2[0.1] == 0.00}, {y1, y2}, {t, 0., tp}, MaxStepSize -> 1/4*tp];

Plot[{Qin[t], Qin[t] - 1/Ithroat Evaluate[y1[t] /. s1]}, {t, 0, 1}, 
 PlotRange -> {{0, 1}, {0, 1.2}}, Frame -> True, Axes -> True, 
 FrameLabel -> {"Time (secs)", "Output flow rates (L/s)"}, 
 GridLines -> Automatic, AspectRatio -> .65, 
 PlotStyle -> {{Red, Thickness[0.01]}, {Blue, Thickness[0.01]}}, 
 BaseStyle -> {FontWeight -> "Bold", FontSize -> 15}]
POSTED BY: Kaushik Mallick
3 Replies
Posted 10 years ago

It is because you require the t in Qint to be numeric? Why? If you eliminate this, you get rid of all the errors. I attached the corrected version using Module to localise all variables, which is safer.

Qin[t_(*?NumericQ*)] := fsum[\[Theta]];

enter image description here

Attachments:
POSTED BY: Erik Mahieu
Posted 10 years ago

It is because you require the t in Qint to be numeric? Why? If you eliminate this, you get rid of all the errors. I attached the corrected version using Module to localise all variables, which is safer.

Qin[t_(*?NumericQ*)] := fsum[\[Theta]];

enter image description here

Attachments:
POSTED BY: Erik Mahieu
Posted 10 years ago

Thank you so much for taking the time to look at my problem and suggesting a solution. I was using ?NumericQ in my function definition because I read it here:

http://mathematica.stackexchange.com/questions/18393/what-are-the-most-common-pitfalls-awaiting-new-users/26037#26037

But your suggestion has solved my problem. Appreciate your help.

POSTED BY: Kaushik Mallick
Reply to this discussion
Community posts can be styled and formatted using the Markdown syntax.
Reply Preview
Attachments
Remove
or Discard

Group Abstract Group Abstract