# NDSolve WorkingPrecision

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 NDSolveIDA, only machine real code is available. Unable to continue with complex values or beyond floating-point exceptions.

Any ideas?

• 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. – Michael E2 Feb 11 '18 at 14:04
• @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. – 407PZ Feb 11 '18 at 14:12
• Yeah, I was testing it out. The piecewise nature of Max gives some trouble when differentiated. I didn't solve it. – Michael E2 Feb 11 '18 at 14:17
• The Max function is actually only to avoid vF[tau] dropping down to 0. There might be another way to do the same thing. – 407PZ Feb 11 '18 at 15:21
• When rationalizing equations, èqs=Rationalize[ ....,0]I don't have any problems with NDSolve, working with MMA 8.0. – Akku14 Feb 11 '18 at 17:19

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.

• 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. – 407PZ Feb 11 '18 at 18:37