I am using Mathematica to solve a system of differential equations and evaluate and plot the solutions on an interval. The solutions have non-zero solutions, but Mathematica underflows in a portion of the expression and as a result provides an incorrect zero result. It seems like there is a way to overcome this but I seem unable to implement it in my code. My code is as follows:


Ac = 1.9793*10^(-6); (*m^2*)
k = 7.1;(*W/m*K*)
P = 0.005;(*m*)
U = \
Cp = 553;(*J/kg*K*)
q = 300/(P*h);(*W/m^2*)
\[Sigma] = 
\[Epsilon] = \
To = 295;(*K*)
Ti = 295;(*K*)
Tl = 2273;(*K*)
h = 0.025;(*m*)
l = 0.05;(*m*)
\[Rho] = 4506;(*kg/m^3*)

eqns = {(Ac*k)*x''[t] - (\[Rho]*Ac*U*Cp)*x'[t] + P*q == 0,
    (Ac*k)*y''[t] - (\[Rho]*Ac*U*Cp)*y'[t] - hr[Ta]*y[t] == 0,
    x[0] == To, y[l] == Tl, x[h] == y[h], x'[h] == y'[h]};
s = DSolve[eqns, {x[t], y[t]}, t][[1]];
Y[t_] = y[t] /. s // FullSimplify
X[t_] = x[t] /. s // FullSimplify

hr[Ta] = 300; NotDone = True; n = 1;
While [NotDone, {Y[t]; Ta[Y] = (1/(l - h))*Integrate[Y[t], {t, h, l}];
   hrn[Ta] = \[Sigma]*\[Epsilon]*P*(Ta[Y] + Ti) (Ta[Y]^2 + Ti^2); 
  err = Abs[(hr[Ta] - hrn[Ta])/hr[Ta]]; NotDone = err > 0.001; 
  hr[Ta] = hrn[Ta];  n++}]

The average value of the function returns a correct result of approx. 1740 but evaluating the function directly returns a zero value and the error "Exp[-877.895] is too small to represent as a normalized machine number; precision may be lost."

Is there a simple way to circumvent this?


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.