Getting message NDSolve::ndnum when solving my initial value problem

Question

I am eagerly learning Mathematica and trying to solve a initial value problem using NDSolve as follows. I'm running into the message NDSolve::ndnum.

ClearSystemCache;
ClearAll["Global`*"];

(*Saturation pressure (Pa) as function of Temperature (K) *)

psatR134a[T_] := 
  Module[{a = -33.7906321345355, b = 0.206544045004122, 
   c = 4.03064679080556, d = 4780.59948872221, p}, p = (a + b*T)^c + d]

(*Saturation Temperature (K) as function of pressure (Pa) *)

tsatR134a[p_] := 
 Module[{a = -732.46, b = 13.733, c = 0.51826, d = 177.01, T}, 
  T = (a + b*Sqrt[p])^c + d]

dqdt[T_, p_, q_] := 
 Module[{A1 = 0.0859, A2 = 0.0566, B1 = 9.14 , B2 = -8.3, 
   Tcri = 374.21(*K*), pcri = 4059.48*^3 (*Pa*), Ea := 7332.69 (*J/
   mol*), R := 8.3144 (*J/K.mol*), qt},
  qt = ((A1*(p/pcri) + A2)*Exp[-Ea/R/T]*(1. + (1. - Exp[-1.5*p/p])) + 
      Exp[B1*(tsatR134a[p]/Tcri) + B2]*(1 - tsatR134a[p]/T))*(qR134a[
       T, p] - q)]

qR134a[T_, p_] := 
 Module[{q0 := 2.2 (*kg/kg*), Ea := 7332.69 (*J/mol*), n := 1.29, 
   R := 8.3144 (*J/K.mol*), qs}, 
  qs = q0*Exp[-(R*T/Ea*Log[psatr134a[T]/p])^n]]

tOutlet[tIn_, tp_, aR_, uH_, mdot_, cp_] := 
 Module[{tAirOut}, tAirOut = tp + (tIn - tp)*Exp[-uH*aR/mdot/cp]]

tOutlet[298.15, 285.5760375366, 6.25, 15, 0.0851, 1006.6]

(* 289.785 *)

ClearAll[tads, tevap, mR134a, qa, var, vars, eqns, ics]

vars := {tads[t], tevap[t], mR134a[t], qa[t]};

eqns := Module[
   {

    mads = 10,
    cpads = 1375,
    cpR134a = 1400,
    mRr134a = 10,
    hgR134a = 420*^3,
    Qst = 325*^3,
    mHXE = 1.958,
    cpE = 910,
    mdotAir = 0.0851,
    cpAir = 1006.6,
    tAirIn = 298.15,
    uAir = 15,
    areaE = 6.25,
    },
   {

    (mads*cpads + qa[t]*cpR134a)*tads'[t] == Qst*mads*qa'[t],
    qa'[t] == dqdt[tads[t], psatR134a[tevap[t]], qa[t]],
    mR134a'[t] == -mads*qa'[t],
    (mHXE*cpE + mR134a[t]*cpR134a)*tevap'[t] == -hgR134a*mads*qa'[t] +
       mdotAir*
       cpAir*(tAirIn - 
         tOutlet[tAirIn, tevap[t], areaE, uAir, mdotAir, cpAir])
    }
   ];

ics := {tads[0] == 303.15, tevap[0] == 303.15, mR134a[0] == 10,
   qa[0] == 0.05}

s = NDSolve[{eqns, ics}, vars, {t, 0, 300}]

During evaluation of this last line, I receive the message,

NDSolve::ndnum: Encountered non-numerical value for a derivative at t == 0.`. >>

I got the solution in Fortran as below.

Your help is greatly appreciated.

What happens if you plug in your ICs and t == 0 into your DE? — Michael E2, Aug 12 '16 at 00:53
Thank you for your suggestion. Here is what I got. NDSolve[{{14450. Derivative[1][tads][t] == 3250000 Derivative[1][qa][t], Derivative[1][qa][t] == 0.00704719 (-0.5 + 2.2 E^(-0.252192 Log[1.29894*10^-6 psatr134a[303.15]]^1.29)), Derivative[1][mR134a][t] == -10 Derivative[1][qa][t], 15781.8 Derivative[1][tevap][t] == -284.939 - 4200000 Derivative[1][qa][t]}, {tads[0] == 303.15, tevap[0] == 303.15, mR134a[0] == 10, qa[0] == 0.5}}, {tads[t], tevap[t], mR134a[t], qa[t]}, {t, 0, 300}] — yaykhel, Aug 12 '16 at 01:05
Just a side note, ClearSystemCache won't work, it should be ClearSystemCache[]. — xzczd, Aug 12 '16 at 05:50

xinxin guo · Accepted Answer · 2016-08-12T05:48:41.623

Function psatr134a is missing. If it is given as a Exp, then a solution can be obtained.

As pointed out by @xzczd in the comments, if psatr134a changed to psatR134a, a solution can be obtained.

ClearSystemCache;

ClearAll["Global`*"];

(*Saturation pressure (Pa) as function of Temperature (K)*)

psatR134a[T_?NumericQ] := 
 Module[{a = -33.7906321345355, b = 0.206544045004122, 
  c = 4.03064679080556, d = 4780.59948872221, p}, p = (a + b*T)^c + d]

(*Saturation Temperature (K) as function of pressure (Pa)*)

tsatR134a[p_?NumericQ] := 
 Module[{a = -732.46, b = 13.733, c = 0.51826, d = 177.01, T}, 
 T = (a + b*Sqrt[p])^c + d]

dqdt[T_?NumericQ, p_?NumericQ, q_?NumericQ] := 
 Module[{A1 = 0.0859, A2 = 0.0566, B1 = 9.14, B2 = -8.3, 
 Tcri = 374.21(*K*), pcri = 4059.48*^3 (*Pa*), Ea = 7332.69 (*J/mol*),
 R = 8.3144 (*J/K.mol*), qt}, 
 qt = ((A1*(p/pcri) + A2)*Exp[-Ea/R/T]*(1. + (1. - Exp[-1.5*p/p])) + 
  Exp[B1*(tsatR134a[p]/Tcri) + B2]*(1 - tsatR134a[p]/T))*(qR134a[
   T, p] - q)]

  qR134a[T_?NumericQ, p_?NumericQ] := 
  Module[{q0 := 2.2 (*kg/kg*), Ea := 7332.69 (*J/mol*), n := 1.29, 
  R := 8.3144 (*J/K.mol*), qs}, 
  qs = q0*Exp[-(R*T/Ea*Log[psatR134a[T]/p])^n]]

 tOutlet[tIn_?NumericQ, tp_?NumericQ, aR_?NumericQ, uH_?NumericQ, 
 mdot_?NumericQ, cp_?NumericQ] := 
 Module[{tAirOut}, tAirOut = tp + (tIn - tp)*Exp[-uH*aR/mdot/cp]]

 tOutlet[298.15, 285.5760375366, 6.25, 15, 0.0851, 1006.6]



 ClearAll[tads, tevap, mR134a, qa, var, vars, eqns, ics]

 vars = {tads[t], tevap[t], mR134a[t], qa[t]};

 {mads = 10, cpads = 1375, cpR134a = 1400, mRr134a = 10, 
 hgR134a = 420*^3, Qst = 325*^3, mHXE = 1.958, cpE = 910, 
 mdotAir = 0.0851, cpAir = 1006.6, tAirIn = 298.15, uAir = 15, 
 areaE = 6.25}; eqns = {(mads*cpads + qa[t]*cpR134a)*tads'[t] == 
 Qst*mads*qa'[t], 
 qa'[t] == dqdt[tads[t], psatR134a[tevap[t]], qa[t]], 
 mR134a'[t] == -mads*qa'[t], (mHXE*cpE + mR134a[t]*cpR134a)*
 tevap'[t] == -hgR134a*mads*qa'[t] + 
 mdotAir*cpAir*(tAirIn - 
   tOutlet[tAirIn, tevap[t], areaE, uAir, mdotAir, cpAir])};

ics = {tads[0] == 303.15, tevap[0] == 303.15, mR134a[0] == 10, 
 qa[0] == 0.05}

 s = NDSolve[{eqns, ics}, vars, {t, 0, 300}]

Hi xinxin guo and xzczd, greatly appreciate your help. It works now. May I know why do you put "?NumericQ" in the function argument? Can get some recommended books to be Mathematica experts like you? Many thanks. — yaykhel, Aug 12 '16 at 06:10
For ?NumericQ: http://mathematica.stackexchange.com/a/26037 -- For general references, look for the "....new to Mathematia" part: is this: http://mathematica.stackexchange.com/questions/18/ — Michael E2, Aug 12 '16 at 10:29

Getting message NDSolve::ndnum when solving my initial value problem

1 Answers1