SummaryWe are trying to solve some models of emission of substances to the air. In these models, emission stops at a certain point.We are using the DotNumerics implementation of radau5We had issue finding correct solutions. If emission stops, the air concentration should drop, but this was not always the case.We were able to reproduce this issue with a simple model for which we could also find the analytic solution. Also, we inspected the radau5 source code and found a test that looks like an optimization. Due to this test many evaluations of the Jacobian are skipped. If we remove this test, the results are a lot better.We have also used the implicit fourth order Runge Kutta by Press et all., Numerical Recipes 2nd edition for our original problems and this method did not show this issue. We didn't find the time to test it on this problem.We wonder:What is this test? Is there some documentation on it?Can we safely remove the test to skip evaluation of the Jacobian?Can we expect radau5 to cope with discontinuities, or should we integrate the two trajectories separately?Detailed informationThe model consists of two compartments. Initially, compartment I contains the substance that is emitted to compartment II. From compartment II, product is removed by ventilation.$$\begin{array}{rclr}\frac{dC_I}{dt}&=&-k_{12}(C_I- C_{II}) \\\frac{dC_{II}}{dt}&=&k_{12}(C_I- C_{II})-k_{02}C_{II} & t < T \\\frac{dC_{II}}{dt}&=&-k_{02}C_{II} & t > T \\C_I (0)&=&1 \\C_{II}(0)&=&0\end{array}$$We found these analytic solutions:$$\begin{array}{rclr}C_{II}(t)&=&C_{II}(T)e^{-k_{02}(t-T)} \\&=&\frac{C_0}{\sqrt{D}}k_{12}(e^{_+ T}-e^{_-T})e^{-k_{02}(t-T)} & t < T \\C_{II}(t)&=&\frac{C_0}{\sqrt{D}} k_{12} (e^{_+t}-e^{_- t}) & t > T \\D&=&4k_{12}^2+k_{02}^2 \\_&=&-k_{12}-\frac{k_{02}}{2}\frac{\sqrt{D}}{2}\end{array}$$We implemented this model in C#.The numeric solution with a tolerance of 0.01 is as follows:jsfiddle for original radau5 and relative tolerance: 0.01In this case, the Jacobian is evaluated only twice.Looking through the radau5 code, we found a test:if (THETA <= THET) goto LABEL20;on line 1663 in .\DotNumerics\ODE\Radau5\radau5.cs, IF (THETA.LE.THET) GOTO 20on line 1073 in http://www.unige.ch/~hairer/prog/stiff/radau5.f.If we apply the fix we get:jsfiddle for original radau5 and relative tolerance: 0.01Reducing the tolerance to 0.001 yields better result, with or without fix.jsfiddle for original radau5 and relative tolerance: 0.001jsfiddle for modified radau5 and relative tolerance: 0.001However, without fix, the Jacobian is still evaluated only three times, so that could be a coincidence.A c# solution that demonstrates this is on GitHub.
Failing integration with the radau5-implementation in DotNumerics, over a discontinuity
numerical analysis;ode;stiffness