1
$\begingroup$

I'm trying to solve differential equations but encountered singularity problem as usual. My idea is to solve is by using higher WorkingPrecision, but somehow it always gives me errors

For example,

eqs = {51.7 ti'[tau] == qH[tau] + qS - 0.7343549132947976 (-tae + ti[tau]), 
  0 == -ti[tau] + tr[tau] - Exp[-2.655/Max[1.*10^(-9), vF[tau]]] (-ti[tau] + tv), 
  0 == qH[tau] - 1.1021142600089968` (-tr[tau] + tv) vF[tau],
  0 == vF[tau] - Max[0, 20 + 20 (20 - ti[tau])]}
vars = {ti[tau], qH[tau], tr[tau], vF[tau]}

The solution of the steady-state equation can be used as initial condition. Without manually setting the working precision, it runs smoothly.

ic = (FindRoot[eqs[[All, 2]] //. {tv -> 28, qS -> 0, tae -> -2.5 Cos[2 Pi/24 tau], tau -> 0}, Transpose[{vars /. tau -> 0, ConstantArray[RandomReal[], Length[vars]]}]] /. Rule -> Equal)
NDSolve[{eqs /. {tv -> 28, qS -> 0, tae -> -2.5 Cos[2 Pi/24 tau]}, ic}, vars, {tau, 0, 24 30}]

When I try to set WorkingPrecision->10, it gives me the error.

With[{eqNew = 
   SetPrecision[{eqs /. {tv -> 28, qS -> 0, 
       tae -> -2.5 Cos[2 Pi/24 tau]}, ic}, 10]}, 
 NDSolve[eqNew, vars, {tau, 0, 24 30}, WorkingPrecision -> 10]]

NDSolve::mconly: For the method NDSolve`IDA, only machine real code is available. Unable to continue with complex values or beyond floating-point exceptions.

Any ideas?

Improve this question
$\endgroup$
8
  • 1
    $\begingroup$ You can only use WorkingPrecision -> MachinePrecision for the IDA method, which is the only way to solve a DAE. Try differentiating your algebraic equations in eqs. -- Side remark: Machine precision is almost 16, but without precision tracking and extra guard digits. I'd recommend setting WorkingPrecision to 16 or higher, if you want higher precision. $\endgroup$
    – Michael E2
    Commented Feb 11, 2018 at 14:04
  • $\begingroup$ @Michael E2 Thanks for reply. This is quite a sad news to me. I struggle with NDSolve all the time. It always gives me error like singularity issues etc. Thought tune up the WorkingPrecision might help. Also appreciate your side remark. $\endgroup$
    – 407PZ
    Commented Feb 11, 2018 at 14:12
  • $\begingroup$ Yeah, I was testing it out. The piecewise nature of Max gives some trouble when differentiated. I didn't solve it. $\endgroup$
    – Michael E2
    Commented Feb 11, 2018 at 14:17
  • $\begingroup$ The Max function is actually only to avoid vF[tau] dropping down to 0. There might be another way to do the same thing. $\endgroup$
    – 407PZ
    Commented Feb 11, 2018 at 15:21
  • $\begingroup$ When rationalizing equations, èqs=Rationalize[ ....,0]I don't have any problems with NDSolve, working with MMA 8.0. $\endgroup$
    – Akku14
    Commented Feb 11, 2018 at 17:19

1 Answer 1

Reset to default
1
$\begingroup$

Reduce the 4 equations to one differential equation with Eliminate.

eqs2 = Rationalize[eqs, 0]

eli1 = Eliminate[eqs2, {qH[tau], tr[tau]}]

eli2 = Rationalize[
       eli1[[1]] /. vF[tau] -> Max[0, 20 + 20 (20 - ti[tau])] /. {tv -> 28,
       qS -> 0, tae -> -2.5 Cos[2 Pi/24 tau]}, 0]

(*   0 == 7060436955/2 Cos[(\[Pi] tau)/12] - 
  59339000000 Max[0, 20 + 20 (20 - ti[tau])] + 
  59339000000 E^(-(531/(
  200 Max[1/1000000000, 20 + 20 (20 - ti[tau])])))
  Max[0, 20 + 20 (20 - ti[tau])] + 1412087391 ti[tau] + 
  2119250000 Max[0, 20 + 20 (20 - ti[tau])] ti[tau] - 
  2119250000 E^(-(531/(
  200 Max[1/1000000000, 20 + 20 (20 - ti[tau])])))
  Max[0, 20 + 20 (20 - ti[tau])] ti[tau] + 
  99413671500 Derivative[1][ti][tau]   *)

ini = Rationalize[ic[[1]], 0]

ndsol = NDSolve[eli2 && ini, ti, {tau, 0, 24 30}]

It workes down to WorkingPrecision->9 with MMA 8.0

ndsol9 = NDSolve[eli2 && ini, ti, {tau, 0, 24 30}, 
         WorkingPrecision -> 9]

You then get the qH[tau], tr[tau], vF[tau] from the algebraic equations.

Improve this answer
$\endgroup$
1
  • $\begingroup$ It's a very nice solution! It would be even nicer, if I won't have to simplify the equations with Eliminate, because what I posted here are only simplified solutions. It can be very exhausted to look up all the equations and try to find out all the simplification possibilities. $\endgroup$
    – 407PZ
    Commented Feb 11, 2018 at 18:37

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.