14
$\begingroup$

I'm trying to make some numerical simulation with NDSolve. I have encountered a few problems. Here is a simplified version of the equations:

Cwirk = 50 25 3; eq= {Q[tau] - 200 (tt[tau] - tae[tau]) == Cwirk tt'[tau], Q[tau] - 100 (28 - tr[tau]) vF[tau] 7/6 == 0, tr[tau] - tt[tau] - (28 - tt[tau]) Exp[-0.9/(7/6 vF[tau] 0.22)] == 0, vF[tau] - Clip[(20 + 20 (20 - tt[tau])), {10^(-10), 100}] == 0, tt[1] == 15, tr[1] == 15, vF[1] == 100, Q[1] == 0}

Using Clip is to limit the vF in a specific range, e.g. between 0 and 100. 10^-10 is to prevent 1/0 error. I don't know if WhenEvent could do the same thing. The tae is an interpolated function. E.g.:

li = Import[http://web.archive.org/web/20180623011154/https://rredc.nrel.gov/solar/old_data/nsrdb/1991-2005/data/tmy3/725958TYA.CSV]; tae = Interpolation[Transpose[{Range[8760], Drop[Drop[li, 1][[All, 32]], 1]}]]

But NDSolve complains about complex values if I try to solve the equations

sol0 = Flatten[NDSolveValue[eq, tt, {tau, 1, 24 365}]] NDSolveValue::mconly: For the method IDA, only machine real code is available. Unable to continue with complex values or beyond floating-point exceptions.

However, if I eliminate some equations with Rule, it runs smoothly.

rs = {Q -> 100 (28 - tr) vF 7/6, tr -> tt[tau] + (28 - tt[tau]) Exp[-0.9/(7/6 vF 0.22)], vF -> Clip[20 + 20 (20 - tt[tau]), {10^(-10), 100}]}; eq = Q - 200 (tt[tau] - tae[tau]) == Cwirk tt'[tau]; sol2 = Flatten[NDSolve[{eq //. rs, tt[1] == 15}, tt, {tau, 1, 24 365}]]

But why?

The real problem is much bigger and it's not very easy to eliminate some equations. Is there any other solution for this problem?

EDIT

Here are part of the data from the file.

li = {{0, 1.5}, {1, -1}, {2, 0}, {3, 0}, {4, 3}, {5, 6}, {6, 10}, {7, 12}, {8, 14}, {9, 16}, {10, 17}, {11, 18}, {12, 20}, {13, 23}, {14, 23}, {15, 21}, {16, 17}, {17, 14}, {18, 12}, {19, 10}, {20, 9}, {21, 8}, {22, 4}, {23, 4}, {24, 1.5}}; tae = Interpolation[li]; sol0 = Flatten[NDSolveValue[eq, tt, {tau, 0, 24}]] NDSolveValue::ivres: NDSolve has computed initial values that give a zero residual for the differential-algebraic system, but some components are different from those specified. If you need them to be satisfied, giving initial conditions for all dependent variables and their derivatives is recommended. NDSolveValue::mconly: For the method IDA, only machine real code is available. Unable to continue with complex values or beyond floating-point exceptions.

I think the Exp[-0.9/(7/6 vF[tau] 0.22)] is the problem. When vF drops to 10^-10, this Exp[] become extremely small near zero and NDSolve has difficulty to process it. Is there any way to avoid this?

EDIT2

@xzczd discovered that Exp[-0.9/(7/6 vF[tau] 0.22)] is not the root of this error message after trying deactivating CatchMachineUnderflow. It appears that NDSolve cannot process the 10^-10 in the eq. Increasing the 10^-10 to 0.1 solves the problem, but doesn't give me the accuracy I want. Any other solutions?

P.S. my intention is to prevent the vF becoming negative, so if there is any way to make Clip[(20 + 20 (20 - tt[tau])), {0, 100}] working is also acceptable.

EDIT3

Here are the optimized initial conditions.

eq = {Q[tau] - 200 (tt[tau] - tae[tau]) == Cwirk tt'[tau], Q[tau] - 100 (28 - tr[tau]) vF[tau] 7/6 == 0, tr[tau] - tt[tau] - (28 - tt[tau]) Exp[-0.9/(7/6 vF[tau] 0.22)] == 0, vF[tau] - Clip[(20 + 20 (20 - tt[tau])), {10^(-10), 100}] == 0} ic = FindRoot[eq[[All, 1]] == 0 /. tau -> 0, Transpose[{{Q[0], tt[0], tr[0], vF[0]}, {3000, 20, 20, 20}}]] /. Rule -> Equal sol = NDSolve[Join[eq, ic], {Q, tt, tr, vF}, {tau, 0, 24}] NDSolve::mconly: For the method NDSolve`IDA, only machine real code is available. Unable to continue with complex values or beyond floating-point exceptions.
Improve this question
$\endgroup$
13
  • 2
    $\begingroup$ Avoid DAE systems, if at all possible. $\endgroup$ Commented Mar 8, 2018 at 2:09
  • 2
    $\begingroup$ You might want to look into the Advanced Numerical Differential Equation Solving Tutorial. There are lots of tweeks to the Methods that may work in your case. $\endgroup$ Commented May 31, 2018 at 12:08
  • 1
    $\begingroup$ I suggest embedding the data for li into your post instead of giving it as an external link, that'll make your post more attractive. I know the original file is too large to post here, but I guess it's enough to reproduce the issue with just a small part of the data? (Actually I have difficulty in downloading the file here… ) $\endgroup$ Commented Jun 1, 2018 at 17:26
  • $\begingroup$ @xzczd thanks for your advice. I added some data into this post. $\endgroup$ Commented Jun 1, 2018 at 18:43
  • 1
    $\begingroup$ Then this is a story telling us to avoid external link in question. Anyway I manage to find a backup, link fixed. $\endgroup$ Commented Mar 29, 2021 at 14:19

5 Answers 5

Reset to default
6
+500
$\begingroup$

For the equations you have given in Case 3, putting a second Clip around the exponential is sufficient to keep the code to machine precision:

Cwirk = 50 25 3; li = {{0, 1.5}, {1, -1}, {2, 0}, {3, 0}, {4, 3}, {5, 6}, {6, 10}, {7, 12}, {8, 14}, {9, 16}, {10, 17}, {11, 18}, {12, 20}, {13, 23}, {14, 23}, {15, 21}, {16, 17}, {17, 14}, {18, 12}, {19, 10}, {20, 9}, {21, 8}, {22, 4}, {23, 4}, {24, 1.5}}; tae = Interpolation[li]; eq = {Q[tau] - 200 (tt[tau] - tae[tau]) == Cwirk tt'[tau], Q[tau] - 100 (28 - tr[tau]) vF[tau] 7/6 == 0, tr[tau] - tt[tau] - (28 - tt[tau]) Clip[Exp[-0.9/(7/6 vF[tau] 0.22)], {10^-10, 100}] == 0, vF[tau] - Clip[(20 + 20 (20 - tt[tau])), {10^(-10), 100}] == 0} ic = FindRoot[eq[[All, 1]] == 0 /. tau -> 0, Transpose[{{Q[0], tt[0], tr[0], vF[0]}, {3000, 20, 20, 20}}]] /. Rule -> Equal sol = NDSolve[Join[eq, ic], {Q, tt, tr, vF}, {tau, 0, 24}]

Using the full data, the solver thinks the system is singular at some point around t=2124, with NDSolve::ndsz (step size is effectively zero):

li = Import["http://rredc.nrel.gov/solar/old_data/nsrdb/1991-2005/data/tmy3/725958TYA.CSV"]; tae = Interpolation[Transpose[{Range[8760], Drop[Drop[li, 1][[All, 32]], 1]}]]; ic = FindRoot[eq[[All, 1]] == 0 /. tau -> 1, Transpose[{{Q[1], tt[1], tr[1], vF[1]}, {3000, 20, 20, 20}}]] /. Rule -> Equal; solfull0 = NDSolve[Join[eq, ic], {Q, tt, tr, vF}, {tau, 1, 8760}, MaxSteps -> ∞]

If we use the StateSpace method, the solution can be found (in 90 seconds for me):

solfull = NDSolve[Join[eq, ic], {Q, tt, tr, vF}, {tau, 1, 8760}, Method -> {"TimeIntegration" -> {"StateSpace"}}, MaxSteps -> ∞] Plot[tt[tau] /. solfull, {tau, 1, 8760}, AxesLabel -> {tau, "tt"}, PlotRange -> All, LabelStyle -> Directive[Bold, Black, Medium], ImageSize -> Large]

enter image description here

This looks the same as that given in bbgodfrey's post where the equations are transformed into a single ODE.

Improve this answer
$\endgroup$
14
  • $\begingroup$ Although when I try the same for the full interpolation function posted, it gets to somewhere around 2130 before failing with NDSolve::ndsz - singularity or stiff system suspected, but I can't see any reason why this is the case. $\endgroup$ Commented Jun 2, 2018 at 20:00
  • $\begingroup$ Thanks for the effort. It's a pity that this trick only works with part of the data, not the full interpolated function. I've struggle with this problem for months. $\endgroup$ Commented Jun 2, 2018 at 20:35
  • 1
    $\begingroup$ There is a weak singularity because of the clipping there. Seems to trick NDSolve (or IDA) into thinking it's stiff. I suppose it doesn't happen with the OP's 2nd, non-DAE system because LSODA is more robust. But my M just crashed, and I hadn't saved my work on this problem. :/ Not that I felt I was close to solving it, so maybe not a big deal. $\endgroup$ Commented Jun 2, 2018 at 21:32
  • $\begingroup$ @MichaelE2 NDSolve has no difficulty to solve non-DAE version. But in the real problem, it is really painful to eliminate all the non-differential equations (maybe even impossible) because the system is so big. Is there any way to make IDA think it is not a stiff? $\endgroup$ Commented Jun 2, 2018 at 22:00
  • $\begingroup$ @407Peezy Sorry, I don't know. $\endgroup$ Commented Jun 2, 2018 at 22:32
7
$\begingroup$

The 10^-10 in Clip is indeed causing trouble, but not in Exp[-0.9/(7/6 vF[tau] 0.22)], because the mconly warning still pops up even after the technique mentioned in this post is used.

So, what's the root of the problem? I'm not sure, nevertheless, the problem can be resolved by changing 10^-10 to something a bit larger:

eq = {Q[tau] - 200 (tt[tau] - tae[tau]) == Cwirk tt'[tau], Q[tau] == 100 (28 - tr[tau]) vF[tau] 7/6, tr[tau] == tt[tau] + (28 - tt[tau]) (10^-10 + Exp[-0.9/(7/6 vF[tau] 0.22)]), vF[tau] - Clip[(20 + 20 (20 - tt[tau])), {0.1(*See here!*), 100}] == 0, tt[1] == 15, tr[1] == 15, vF[1] == 100, Q[1] == 0}; sol0 = NDSolveValue[eq, tt, {tau, 0, 24}, MaxSteps -> Infinity] {sol0, sol2}[24] // Through (* {20.2045, 20.2045} *)

The result at tau == 24 is the same as that of sol2, so I think this solution is acceptable?

Improve this answer
$\endgroup$
3
  • $\begingroup$ Thanks for your answer. I've also considered your suggestion. However, the error is bigger than my acceptable value. Is there any way to keep the 10^-10 and still get it working? $\endgroup$ Commented Jun 1, 2018 at 19:15
  • 4
    $\begingroup$ @407Peezy Sadly I don't know. Actually similar problems have appeared in this site for many times and I never see them being solved in a direct way i.e. by adjusting options, as you've mentioned. I myself tend to take it as the current limitation of DAE solver in Mathematica. I sincerely hope I'm wrong. $\endgroup$ Commented Jun 1, 2018 at 19:21
  • $\begingroup$ That would be really quite a pity if Mathematica DAE solve is not able to process this kind of equation. $\endgroup$ Commented Jun 1, 2018 at 19:27
7
$\begingroup$

The block of code in the question following "However, if I eliminate some equations with Rule, it runs smoothly" works, because the transformation changes the differential-algebraic system to a purely differential system, for which NDSolve uses more robust solution methods. The same can be accomplished more generally by expressing the variables appearing only algebraically as functions of tt. Here, the functions are

Qf[tt_?NumericQ] := 100 (28 - trf[tt]) vFf[tt] 7/6; trf[tt_?NumericQ] := tt + (28 - tt) Exp[-(270/77)/vFf[tt]]; vFf[tt_?NumericQ] := Clip[(20 + 20 (20 - tt)), {10^-10, 100}];

and the system, now a single ODE, is solved by

Cwirk = 50 25 3; eq = {Qf[tt[tau]] - 200 (tt[tau] - tae[tau]) == Cwirk tt'[tau], tt[1] == 15}; sol1 = Flatten[NDSolveValue[eq, {tt}, {tau, 0, 24}]]; ul = sol1[[1]]["Domain"][[1, 2]]; Plot[sol1[[1]][tau], {tau, 0, ul}, AxesLabel -> {tau, tt}, LabelStyle -> Directive[Bold, Black, Medium], ImageSize -> Large]

enter image description here

In general, the variables appearing only algebraically can be converted into functions, although some care must be taken, if they are not single-valued.

As an aside, the original system of equations in the question, when rewritten as

Cwirk = 50 25 3; eq = {Q[tau] - 200 (tt[tau] - tae[tau]) == Cwirk tt'[tau], vF[tau] == Clip[(20 + 20 (20 - tt[tau])), {3 10^-3, 100}] , tr[tau] == tt[tau] + (28 - tt[tau]) Exp[-0.9/(7/6 vF[tau] 0.22)], Q[tau] == 100 (28 - tr[tau]) vF[tau] 7/6, tt[1] == 15, tr[1] == 15, vF[1] == 100, Q[1] == 0}; sol0 = Flatten[NDSolveValue[eq, {tt, tr, vF, Q}, {tau, 0, 24}]];

also runs to completion, because the lower limit of Clip now is set to 3 10^-3, which I would think would be good enough for most purposes. (The exponential in the question can assume a minimum value as small as about 10^-305.) Incidentally, the warning messages, NDSolveValue::ivres returned by this latter block of code indicate that the NDSolve has corrected the initial conditions for consistency:

Table[sol0[[i]][1], {i, 4}] (* {15., 27.5521, 100., 5226.02} *)

Using these values instead of those in the question eliminates the warning messages.

Addendum: Result for imported li

The first block of code in this answer also can handle the imported version of li, despite the fact that the resulting tae is very noisy.

li = Import["http://rredc.nrel.gov/solar/old_data/nsrdb/1991-2005/data/tmy3/725958TYA.CSV"]; tae = Interpolation[Transpose[{Range[8760], Drop[Drop[li, 1][[All, 32]], 1]}]]; Plot[tae[tau], {tau, 1, 24 365}, AxesLabel -> {tau, "tae"}, PlotRange -> All, LabelStyle -> Directive[Bold, Black, Medium], ImageSize -> Large]

enter image description here

Then, the first block of code, integrated over {tau, 1, 24 365}, yields

enter image description here

Incidentally, the second block of code with the lower bound of Clip set to 5 10^-2 reproduces the final plot above.

Improve this answer
$\endgroup$
8
  • $\begingroup$ Thanks for your answer. However, like already mentioned, the real problem cannot be easily rewritten into pure ODE system. Thus, it would be better, if it doesn't require to transform the DAE system into ODE system. $\endgroup$ Commented Jun 3, 2018 at 0:01
  • $\begingroup$ And the to 2nd part. Although it runs to the end with 3 10^-3 in the simplified version, it still fails to run to completion for the full interpolated function tae[tau]. mconly still appears. $\endgroup$ Commented Jun 3, 2018 at 0:03
  • $\begingroup$ However, I have to admit, it seems that the ODE solver is really more stable than the DAE solver. $\endgroup$ Commented Jun 3, 2018 at 0:05
  • $\begingroup$ @407Peezy Don't be too quick to dismiss using functions to eliminate algebraic equations. How many differential equations and how many algebraic equations are included in your complete problem? $\endgroup$ Commented Jun 3, 2018 at 4:27
  • $\begingroup$ There are 12 differential equations and 12 algebraic ones $\endgroup$ Commented Jun 3, 2018 at 10:48
2
$\begingroup$

I would like to share some new findings to this post after nearly 2 years.

As described in the question, running the original code gives me the NDSolve::mconly error.

li = Import["http://rredc.nrel.gov/solar/old_data/nsrdb/1991-2005/data/tmy3/725958TYA.CSV"]; tae = Interpolation[Transpose[{Range[8760], Drop[Drop[li, 1][[All, 32]], 1]}]] Cwirk = 50 25 3; eq0 = {Q[tau] - 200 (tt[tau] - tae[tau]) == Cwirk tt'[tau], Q[tau] - 100 (28 - tr[tau]) vF[tau] 7/6 == 0, tr[tau] - tt[tau] - (28 - tt[tau]) Exp[-0.9/(7/6 (vF[tau]) 0.22)] == 0, vF[tau] - Clip[20 + 20 (20 - tt[tau]), {$MachineEpsilon, 100}] == 0} ic0 = FindRoot[eq0[[All, 1]] == 0 /. tau -> 1, Transpose[{{Q[1], tt[1], tr[1], vF[1]}, {3000, 20, 20, 20}}], MaxIterations -> 10000] /. Rule -> Equal sol0 = Monitor[NDSolve[Rationalize[Join[eq0, ic0]], {Q, tt, tr, vF}, {tau, 1, 24 365}, EvaluationMonitor :> (monitor = Row[{"t = ", tau}])], monitor][[1]] (* NDSolve::mconly error *)

But if I modify the third equation tr[tau] - tt[tau] - (28 - tt[tau]) Exp[-0.9/(7/6 (vF[tau]) 0.22)] == 0 into tr[tau] - tt[tau] - (28 - tt[tau]) Exp[-0.9/(7/6 Max[$MachineEpsilon, vF[tau]] 0.22)] == 0, NDSolve returns a different error message NDSolve::ndsz. This is quite interesting. My guess is that NDSolve thought it might be a good idea to try vF[tau] == 0 in the equation at some step which might result in NDSolve::mconly error for IDA solver.

The next thing I tried is to eliminate the second equation Q[tau] - 100 (28 - tr[tau]) vF[tau] 7/6 == 0 and surprisingly it works.

eq2 = {100 (28 - tr[tau]) vF[tau] 7/6 - 200 (tt[tau] - tae[tau]) == Cwirk tt'[tau], tr[tau] - tt[tau] - (28 - tt[tau]) Exp[-0.9/(7/6 Max[vF[tau], $MachineEpsilon] 0.22)] == 0, vF[tau] - Clip[20 + 20 (20 - tt[tau]), {$MachineEpsilon, 100}] == 0} ic2 = FindRoot[eq2[[All, 1]] == 0 /. tau -> 1, Transpose[{{tt[1], tr[1], vF[1]}, {20, 20, 20}}], MaxIterations -> 10000] /. Rule -> Equal AbsoluteTiming[sol2 = Monitor[NDSolve[Rationalize[Join[eq2, ic2]], {tt, tr, vF}, {tau, 1, 24 365}, EvaluationMonitor :> (monitor = Row[{"t = ", tau}])], monitor][[1]];] (* 64.5397 *)

This leads me to take a deeper look into the second equation. Actually for the equation solving, the second equation is not necessary. We can solve the tt and calculate the Q afterwards. After some brief investigation, it turns out that the second equation is quite stiff. If we insert the forth equation and the third equation in to the second, we will have a function $Q(tt(\tau))$. If we plot the function $Q(tt(\tau))$, we can easily see where the problem is. 

rl1 = {Q[tau] -> 100 (28 - tr[tau]) vF[tau] 7/6, tr[tau] -> tt[tau] + (28 - tt[tau]) Exp[-0.9/(7/6 Max[$MachineEpsilon, vF[tau]] 0.22)], vF[tau] -> Clip[20 + 20 (20 - tt[tau]), {$MachineEpsilon, 100}]} Plot[(rl1[[1, 2]] //. Drop[rl1, 1]) /. {tt[tau] -> x}, {x, 15, 25}, AxesLabel -> {"tt", "Q"}]

Plot

The curve becomes very steep when the $tt$ increases to 21 and after 21 the curve becomes flat. If we extract the values for the last steps in sol0, we will see that the last values for $tt$ is around 21.

After all these back and forth and trying for really couple years, I come to a conclusion concerning using NDSolve, which I want to tell myself repeatedly and share with everyone:

Always check the equation first before using NDSolve, even they have physical meaning.

Improve this answer
$\endgroup$
2
  • $\begingroup$ Er… then why do you include the second equation at the beginning? $\endgroup$ Commented Feb 12, 2020 at 11:14
  • $\begingroup$ @xzczd Because this variable is something I would like to have at first. $\endgroup$ Commented Feb 12, 2020 at 12:32
1
$\begingroup$

I found out that setting the WorkingPrecision to 8 simply solves this problem.

li = Import["http://rredc.nrel.gov/solar/old_data/nsrdb/1991-2005/data/tmy3/725958TYA.CSV"]; tae = Interpolation[Transpose[{Range[8760], Drop[Drop[li, 1][[All, 32]], 1]}]] Cwirk = 50 25 3; eq = {Q[tau] - 200 (tt[tau] - tae[tau]) == Cwirk tt'[tau], Q[tau] - 100 (28 - tr[tau]) vF[tau] 7/6 == 0, tr[tau] - tt[tau] - (28 - tt[tau]) Exp[-0.9/(7/6 vF[tau] 0.22)] == 0, vF[tau] - Clip[(20 + 20 (20 - tt[tau])), {10^(-10), 100}] == 0} ic = FindRoot[eq[[All, 1]] == 0 /. tau -> 1, Transpose[{{Q[1], tt[1], tr[1], vF[1]}, {3000, 20, 20, 20}}]] /. Rule -> Equal workingPrecision = 8; AbsoluteTiming[ sol = NDSolve[SetPrecision[Join[eq, ic], workingPrecision], {Q, tt, tr, vF}, {tau, 1, 24 365}, WorkingPrecision -> workingPrecision, MaxSteps -> Infinity];] (* {31.2564, Null} *)

But somehow the Rule method doesn't work anymore on my MMA 12.0 mac.

rs = {Q -> 100 (28 - tr) vF 7/6, tr -> tt[tau] + (28 - tt[tau]) Exp[-0.9/(7/6 vF 0.22)], vF -> Clip[20 + 20 (20 - tt[tau]), {10^(-10), 100}]}; eq2 = Q - 200 (tt[tau] - tae[tau]) == tt'[tau]; AbsoluteTiming[ sol2 = Flatten[ NDSolveValue[{eq2 //. rs, tt[1] == 18.887256814210946}, tt, {tau, 1, 24 365}, MaxSteps -> Infinity]]] NDSolveValue::ndcf: Repeated convergence test failure at tau == 3497.9989465540943`; unable to continue

EDIT

Here is another solution for this case. I always thought that NDSolve does not particularly like discontinuous functions such as Clip, Min, Max. It is possible to replace them with some events. Below a discrete variable mode is used to indicate if the function is below zero or above 100. If below zero, then the mode will be set to -1. If above 100, it will be set to 1. Otherwise, it is 0. Using the mode variable, Clip can be avoided.

Cwirk = 50 25 3; li = Import["http://rredc.nrel.gov/solar/old_data/nsrdb/1991-2005/data/tmy3/725958TYA.CSV"]; tae = Interpolation[ Transpose[{Range[8760], Drop[Drop[li, 1][[All, 32]], 1]}]] eq1 = { Q[tau] - 200 (tt[tau] - tae[tau]) == Cwirk tt'[tau], Q[tau] - 100 (28 - tr[tau]) vF[tau] 7/6 == 0, tr[tau] - tt[tau] - (28 - tt[tau]) Exp[-0.9/7/6 vF[tau] 0.22] == 0, vF[tau] - DiscreteDelta[mode[tau] - 1] 100 - DiscreteDelta[mode[tau]] (20 + 20 (20 - tt[tau])) - DiscreteDelta[mode[tau] + 1] $MachineEpsilon == 0} ic1 = Join[ FindRoot[eq1[[All, 1]] == 0 /. {tau -> 1, mode[_] -> 0}, Transpose[{{Q[1], tt[1], tr[1], vF[1]}, {3000, 20, 20, 20}}]], {mode[1] -> 0}] /. Rule -> Equal evts1 = {WhenEvent[(20 + 20 (20 - tt[tau])) < $MachineEpsilon, mode[tau] -> -1], WhenEvent[(20 + 20 (20 - tt[tau])) > 100, mode[tau] -> 1], WhenEvent[(20 + 20 (20 - tt[tau])) < 100, mode[tau] -> 0], WhenEvent[(20 + 20 (20 - tt[tau])) > $MachineEpsilon, mode[tau] -> 0]} AbsoluteTiming[ sol1 = NDSolve[ Join[eq1, evts1, ic1], {Q, tt, tr, vF, mode}, {tau, 1, 24 365}, DiscreteVariables -> mode][[1]];] (* {73.7385, Null} *)

It's worth mentioning that this method is 2x slower than simply setting WorkingPrecision to 8.

Improve this answer
$\endgroup$
3
  • $\begingroup$ Seems that "StateSpace" method is triggered by WorkingPrecision option. ({state}=NDSolve`ProcessEquations[…]; state["NumericalFunction"]["FunctionExpression"] shows that NDSolve is still dealing with a DAE system, StateSpace seems to be the only available DAE solver for arbitrary precision in this case, and adding Method -> {"TimeIntegration" -> {"StateSpace"}} to your code doesn't seem to change the timing. ) It's a bit surprising that adding WorkingPrecision->8 speeds up the code. $\endgroup$ Commented Nov 15, 2019 at 12:56
  • $\begingroup$ Another somewhat surprising thing is, adding SolveDelayed->True to the sample fails in v12 resolves the problem. $\endgroup$ Commented Nov 15, 2019 at 13:03
  • $\begingroup$ @xzczd I guess that the speeding up with WorkingPrecision->8 is probably because the default AccuracyGoal and PrecisionGoal is half of the WorkingPrecision according to the NDSolve documentation. $\endgroup$ Commented Nov 16, 2019 at 16:19

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.