# Getting message NDSolve::ndnum when solving my initial value problem [closed]

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[
{

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,
},
{

(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.

• What happens if you plug in your ICs and t == 0 into your DE? Commented Aug 12, 2016 at 0: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}] Commented Aug 12, 2016 at 1:05
• Just a side note, ClearSystemCache won't work, it should be ClearSystemCache[]. Commented Aug 12, 2016 at 5:50

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,
mR134a'[t] == -mads*qa'[t], (mHXE*cpE + mR134a[t]*cpR134a)*

• I guess it should be psatR134a :D Commented Aug 12, 2016 at 5:43
• For ?NumericQ`: mathematica.stackexchange.com/a/26037 -- For general references, look for the "....new to Mathematia" part: is this: mathematica.stackexchange.com/questions/18 Commented Aug 12, 2016 at 10:29