4
$\begingroup$

I have a very large system of ODEs (up to several thousands ODEs and up to several hundreds of thousands of various coefficients, which power these ODEs). The ODEs and coefficients are generated based on some model parameters and statistical distributions of coefficients. As the number of ODEs and coefficients increase, the probability that Mathematica will fail with "NDSolve::icfail: Unable to find initial conditions that satisfy the residual function within specified tolerances. Try giving initial conditions for both values and derivatives of the functions." substantially increases. Several consecutive runs with the same model parameters, which result in the same number of equations but different values of coefficients may work or may fail. I wonder if anyone has any ideas what to do with that.

Here is and example of a failing ODE system: https://drive.google.com/file/d/0B92m3iuLY2XPZ1ZCUVBsUllSWkE/view?usp=sharing . It is too big to post here. The example fails in 11.2 but works in 10.0!!! However, 10.0 still fails for larger systems.

Thanks a lot.

UPDATE

Ok. I looked at what was going on under the hood and the issue seems to be with Method -> {"EquationSimplification" -> "Residual"} . Once this is removed, then everything works as a clock. The problem is that I personally set this method a while ago because for some values of parameters the resulting system is DAEs rather than pure ODEs and that method was the only, which worked...

Improve this question
$\endgroup$
8
  • $\begingroup$ If you don't have good starting values for your parameters/coefficients to begin with, then you're pretty much doomed, no? $\endgroup$
    – J. M.'s eventual burnout ♦
    Oct 15, 2017 at 3:00
  • $\begingroup$ It is a first order ODE system with just initial values and no constraints. Give it anything as an input and it must solve it. But it does not! $\endgroup$
    – Konstantin Konstantinov
    Oct 15, 2017 at 3:05
  • $\begingroup$ Haven't looked, but if the system is linear and first-order, you should be able to use MatrixExp[]. If it's nonlinear... good luck, dude. $\endgroup$
    – J. M.'s eventual burnout ♦
    Oct 15, 2017 at 3:07
  • $\begingroup$ Of course, it is non linear and even if it were linear MatrixExp is doomed for very large matrices. The question is how to tweak various parameters spread out inside 300+ pages of documentation for NDSolve, so that to coerce it to work. $\endgroup$
    – Konstantin Konstantinov
    Oct 15, 2017 at 3:10
  • $\begingroup$ My point was that for nonlinear equations, you really don't have a guarantee of "niceness", in that if it works, consider yourself lucky (cf. the Lorenz equations). $\endgroup$
    – J. M.'s eventual burnout ♦
    Oct 15, 2017 at 4:37

3 Answers 3

Reset to default
5
$\begingroup$

There exist a number of posts related to this issue in this site, for example:

What's behind Method -> {"EquationSimplification" -> "Residual"}

Why does NDSolve need to solve for the derivatives if the equations are already explicitly solved?

Introduction to Vectors and NDSolve with System of Equations

Plot with different colors the results of an NDSolve with vector solution

Numerically solve large ODE system

System of ODEs specified by large coefficient matrices

There should be more. In short, large symbolic ODE system is burdensome, it's better to vectorize the system or make small black-box function if possible. In your case, the system can be transformed to the following form:

$$U'(t)=m_0+m_1.U(t)+m_2.U(t).U(t)$$

According to your statement here, I guess you already have $m_0$ etc. at hand, so it should be easy for you to build such a equation, but here I'll rebuild this equation from the code given by you with the help of powerful CoefficientArrays:

{ode, ic} =   Uncompress@"1:eJy9nXdgFdW2xpNQpEkISFWKSFBAkRhCUUxO6CEkAYL0GkhCIJQUEFApGkpURKTKNaBE6U0UQcSEoIBIFQiEegEBBURB6hWEp55598zMW+s7a/aZ8+4fXJU9vzOz91prt7W//Xjc8E4Jfj4+PmmF//qjfVJaekJdw7/983dF/vqjVcrI/skJvn/9k+uPtBJ//dEyPjVpVP/0pFHxSf/8t7+fSO2RVuyv/0uvFRQU1LRpcKqP9j9VXHcDLqGBr/aKHZJHprmYnZOGxqelli/6waoDp0PCElyP9zc+bnpg65HjyVc21Qpz/fcOw1+OT2WfT/IzfAn7q+HGXy1kfKD5kWEBs0Jag9dkSZbfP9z9+/v5PujQsf9L+teJg79KPNAcPhCbtbZ8bP9h+gcGWH2gBXxg0derVgaOCdQ/MBC2wcBKF0ftm1MGvJKU1OfV2h9lXy4D3pUluf+IllZ/mv0ITCKqg/2IlvAjCkXu/3Ti9MF6UrzVB1rBB/pUXu2IfyNe/0AC/Lj+rX+Z8HVKgEPib5j0Ysb0n8qcnQxcXUq6UzWoYutLZ4DTSUkrouu/3i490CExAkz6cn9uwLaKD4ATSUl7boXXaHf3OQdvBFLSJkeXiBMbGwLrYEkms7k4ZlvDCzvLOvhYmeA+VhK219qqxbC2h0mEFbO2h0mE7TVXIxG2x4YNTCJsr6UaibAY1vYwibBi1vZae2h7rd3b3vLRVT//tUQn/YclwveflZPjt7PVJpHtYdLdXZVjNp7cKYp7mHT6t4CP5ywOdvBuy5JM1fHL7KeTShy8BsY+iUp12sZqTbB1iklE67B1iklEnbZWI+VEXX193yu+Dr4RWJLl1mnjvnVGr7rTLCA6Vl9Fg+D7F/i1bLgrpLOodTDpUnzazsMvdhW1DiblLCl1vkOF46JeFZNGfx1eteX+E6IYiUk+my9Eb3y7isgLManDe3X63Zr5AwhKUlKT4jMm9fPbB1yRJbk3m7ZWG5s1G0wiDJA1G0wiGps1G0wiDJA1G0wizIYNNJhENDZrNphEGCBrNm2h2eS3P/xE3T9C9aQk+NNDp10bs7DszVDJGB6TLq8fd8i/y209iR2RYVKho/GZHeeMFU0JMalOi9KZw6NfEU0JMelM0evTlox5WBRtkqw2UYTVimWbCJOIxmabCJOIimWbCJOIxmabCJOIJmI9OwI20dCdfYvvfr6t/p0Gw58+2fbixyP3VAmXhExM+mDBmF/f+7JcuKSJMGnbE6NnZfSpFC6pWExakrV10+xVl/UVywZfTKpRJaZYfqka+ndiZyOY9OZDJb+ZuaymyB8xKWRRtYa+u6qFS8wGk3q27ZyTVuoxPYkN45g0Y1zdF3/sWUhPaqtG2lRsSv1KGyuE87bPktw7RTurTsGOSDCJcAo2AmIS4RRs3MIkwinYsQ0mEU7BTvYxiTBl1ikwiXAv1ikwiXAKdpSESYRTDFIjEU6RpEaqeOjtDscDGofzXsSSTO5Vr83auYsmGsbvQ+BPv7s29+MKb50TDQswaUJgs2drnzsvGhZgUpnQoKmBx2aK1q0wqed7P94auuJ90ZopJg1cGFXx+dWy1RpMapBfve+tm1XAwqmUVD02dv2gz2qH89bBktybTaTVxmbNBpMIA2TNBpOIxmbNBpMIA2TNBpOIxmbNBpMIA2TNBpMIs2GjTSQ0m67TvvmpyAnD8kIy/OnE8M83308JEnkRJmWcrDS3RO2HwIKYlLSi0MK6G+o/BNZIpKT3N6+tXeNaoMgfk61WbHurFctaByYR1cFWLCYRTcRWLCYRFctabHtYscEH21dc92w1/TsNhT/9+Y1rbb9JOibancSkpNCO07aMzRftNWHS0kGPJi+eW0nkRZjUbUePZt3HHBU1NiZ90zPo5em9DogaG5M292q++8PC/iIvGmq1saOsNhHb2JhEmA3b2JhENDbr2ZhENBHb2JhEmA3b2JhENDbr2VGwsacvCZtQf1V7/TsNgz/9U+sbKztfOSha6sakb7vezlj47SHRUjcmPbHz+RHZy/8QLXVj0u5Kb8XU/KiwaLUFk25lrO9ec1GGaFiFSZ/vev6ZbgsyRMMqTFqUX/a1gvdkERCTds59dOoXfYeDRTopaXDe3iD/AFnv78YK1vw8bsPmq/oh8RA10s1idar3/fl3PSlSSHLvXtFWnYJ1L0wiHJV1L0winIJ1L0wiHJV1L0winIJ1L0wiHJV1L0wi3IvtczCJcArWvTCJcFTWvTCJcArWvdxYwf91VNa9oqF7+dVzjJxfUFX/dcPhT2euq3aqa8Et0cQSk6K2fji3158lRMMCTKr4cIzjzOKSomEBJs0PWxAfeqSsKJYOt1qxMVYrlrV9TCKqg61YTCKaiK1YTCIqlvXHGFixiTOO/9G5fKL+nUbAn64zNqSewydUNLjGpLTqk+dPGvOCaHCNSeUufxvgf6OoaByBSdse/+jUD/VqiSwWk1p89GaNuuUM2VBseMKkwH2Tfxl7bCNY1WNJ7hu7g9UmYhsbkwizYRsbk4jGZj0bk4jGZr0Ik4gmYhsbkwizYRu7A2zslQ98Dh2vkKwnpcCf/uii7+oSla/ru0F2qIdJfosi79VrXlg01MOk2YVaBoRWla3XY9JjZwraF4nsCZJhpKTLPyb+/NbG3iBnQ0ratOTL6YU+fEIUbTDJP31DMb9K48HAQ0qqW3Oq7zu+E0FPy5LcG2BHq2bDGiAmEabMGiAmEQbIRhtMIsyGNUBMIkyZNUBMIgyQjYCYRJgNa4CYRJgya4AdoQHui52T1S8yVU9KhT99rvCdjWd7rRNN4jBpe4XBl8Zc/1o0icOk4IBXL6Qsl62RYFKfTTs+8h95IlRiypg04bOQjFlRcaJkR0x6pGbfj+df7iVKdsSkJxouHZe02LCVyQ7TMWlNodORH3/5o57EDtMxqaf/71PeHVFKFN8x6X6Xm4sPVDcsaQ9VIy3acCrrcL6BFKVGmns7s+KFr7aCAYaUdG9y1JWhId+BnoIluXf5TlYdlXV5TCKCB+vymEQ4KuvymEQED9blMYlwVNblMYkIHqzLYxLhqKzLYxIRPFiXxyTC5dkeFZMIR2VdHpOI4MG6PCYRjsq6PCYRwYN1+U7Q5TNrP95ldJBhjpoGf/qLK+NDuo27Klpzw6RLx/Yv6dTzSVF8x6TUizdfOpdiMMBhaqTjgx406/XgjJ4UrUbaUim47sv77+lJ7JAOk3pkXF5WsPdPPYkd0qVZbexYq43NRkBMIhqb9WxMIpqIbWxMIsyGbWxMIpqIbWxMIsyGbexY2Ni7vuw8dcGkcH3bpcOfrjbi+d4lu9YXjZUxqUG32wF7R9UUdXeY9HX4sCk7hlQVdXeYlN9mT0Sda+VE0QaTahVkns374KKoa8GkcWen7sne95Ooa8GkWZ+Ftro3/oLIKTBp8r6SP8zMuSZyCkya4R+2NDRkg6i7w6Q1vV8PafrJp6LuLt2qU3S26hRsBMQkwpRZp8Akwr1Yp8AkwinYqIxJhCmzToFJhHuxToFJhCmzToFJhHuxToFJhCmzToFJhHuxTtEZOsXuCtcnT1uyJ1f3QHf408uXDHj17LJ6W3QPsEuNlknsUiMmfZcffaRP5Q36j2A3WyyT2M0WTPrG/3xgt7KN9F/HplpYJrGpFph0x6/6mhOjAvUkdspumcRO2TGp0aO16vg+01df4+zIwzKJDdeYdOuNXosCv3xR/3VsuLZMYsM1Js0PPLknbPELehI7ZbdMYqfsmFRpV+T5ChPP6mucXZi3TGIX5t183a9Pli332nk9ie3WLJPYUSMmJcx/N+pfuxvra5zd5LRMYjc5ManEvR1xIQ066r8u2TZSezVS6yLBby7eFqH/OnYgYZnEDiQwaXvZchvqr4vUk9iBhGUSO5BwU+NHVwZtCbyqr3F2ymmZxE45MelUhbhlwR9s15PYwY1lEju4waT2EcFhEXHh+hpPtY3USY30c1Jxx4e/x+q/Ls02Uqwa6cBD/w79ul8l/del20bqLCT9d2BqTRTUNGQs8/dfuBEF1deYRJ6zkf4Bz+Q5c63Jc5pe53zi50cuDfEVjZixnpkmcwNqgh1mmEiazI2IhOVbnAfdS+RJSPi0/nPZXUtM/PaQ6J1w/vj8Q2PejrhyUETCCa7L+s/zbVItKUzsSm+vvdp45WyZXCAm+ZyrGbrwkUmipCJM2jY2f+GGt1eIMjoxqfKDB2MP3FgiyujEpDF1HhscfayOg29TKenIKxkhA/c+5eDb1O7g1dxC8NKkdsHU1DxH+0fb8QfRXBZHC006J0dCEukkiUhYpsYpZlEoT0LC0UI7yi16J5GWgIiEjyhrB/FEbSc6vyki4QNm1TtldnvqpTh5BDt9c9Cat1sYHlCcHh77fsDY0uP6ifK2MKlK5tXLsX0MgtuKU7HhCV0zc/5sKco3xaSsy738G3ap6+DbVEqa9MYrsUlBDRx8m9odwUxaSKV8QATTpLz1ptgCfo4mag0ecP2FSONbRMLy2E6d2/JbJCQcVTXxL9E7iZTeRCQcVZ3CPsXyJCSsc3SvwdXizR3j5dGiX+6xYe+sGiA6+olJ85rVavZTOwNJccH0fukVl8tnHbAhgq29Mm/prOoHVCKYomeaPrkY8kxNbj5X7AVOjebfwANSL3AKQB0XrW7jflw73CgiiU4hi74Oj/ozd3XIebB6rNwLxlb+I6TVoSAbRrP5Pj023yv55P/jaLaVBYvTpP71VdwKfs5voweHBt+tKdpxwBbnFAorkychYYvTjq6Cj5BanHZ0VUTCFjdhScLk/Owhcos7n5eeGZhxOMxzi2syZ4FPxFQDycsWZ2roGsjiNL1+sD1k+hxNr1+0+INJJ2vVqvvuy3as29g3p7NvHGNfX6Dp9etJinr9K9Ir+/VXixYms9H0+rdY0+s3z4T+kWV/Ygv/q66/wLLsTt3lEiISXm5z6rz+KdouxYNGp15i6TwJSSJHWlb0dVhRUNMdE30dntRPixl0yLFa1nZ4Uu/UqHpcRMJ6UE5NnGIiEtafcSodyEhYC8B5bvyCqMbxGW3nOdVbIhI+E+o8gyZrO3y26vSvFR8asVxWTzh/25lXe1r0dThf1JmfVk30TjgLqNO2T9ZOfsfCZC39g1IlLyUYunrFxemjM64FRe4xkBQXpxMvzUmtnvqm6EAKJo0/W8qvIHia6EAKJo1uM3BSs8dG27DclHZ+36ZjVQwk6XKT4uDK1BE94gMGV9rlKXqjFl2eIhpciS5PEW3RSC5PsaOXdl6eck3/Tp5dngL2DanLU0yv4xQPLpXH/6q0g9YOM4F2dv2F6CSkiIQDqpZcKyKJctdFJBxQu01bsHDkxjR5QO17a/Xy/PuVRRfTYVL28TXzhiQabuZQDKhryld6dtg4w7q0Yhpc+xPbe2f51XNIKtZNAlT3rzqvfmchOJksJeX5ZHc9cCcbHHiUkp4+k9b07FcTwGkLKanLkFXb7/gYdmq9nD5hmiPAIK9dpqN3D9HFNfoHPLu4RjTHlFxcY8cESru4RhTkRRfXiOIOJjnF3a+LcjtxpNdkcEFju/5CpFQsIuFpj3ZsNUdCEp15F5HwIcqkahvXHaqSIo/0q5c2bTsnt7ZIg89NAtj2kFeXxxmuXFTc58gqGHEmu41B80sx0h893Pdk9+b/EZ1PwaQR358LGhpsOK6jmKCa/cyKmX8+bSApJqhu3l2w9+zOdjb0PnuPrGtSNTFKpfdRjPSmNZy//4WN9NrlSiAd22uLSc5LJYrm8Q9II5gmGyPKKRcpIIlIWFilzu7dl4O7JcujRfwq3xmr35MdpHMztpiwNW/jJ+tFZ4bssTjTCl1FZHHatTgg2d59tqBH1+IYMok8uhbHsHvu0bU4hl1Rj67FMex2eXQtjmGZ16NrcexY5iVmth5di2MYOXl0LY4hgnl0Lc6TepL0WhwTyXnlRYU8/gHXX+DlcOdVBOVFJDwudEqfP9giIeEct4hVhRv4B8jeyU2+7z/yoH4iElZwdUokPiIiYRVLp9RYOREJ92NOiZV7ohrH42enpENF0TthhQvnMeRCIhI+3x367yIVs2olgr5VsR8zrbPCHDftIiF9Pya6tEe05ym6tEeUQ4BJL/z2XI9gW7a8iLgjvWrHRNJEq3P4B6Q+rolWi0g412JixsSpcectzP1SQ1q88k7T+zZsUTx+cGHm/b63bdiiiFraeOlrJwy5tYr5HwlVVyZdP2rQafRy/odpGlcYeaZ2pRA4pui17Vqiz1S8vufjQQ9yvz3Y2f4YF2UhxmnXyuhrUnKFix3b1doVLvqf9uwKF9GqlOQKF8M4RnqFi3kz3ilcBipW2j9rwmUiEl7fevi7bw5mjrKwNbypS9G908sZ8g4Vzwh1zx67+ZqvIZtd8YzQ+l2J8eFzDWeEFNdtHpw80fi3PMMZIS+v25hCqD/yTO2SEn2z46QNH5/3Pl849pS+L1S8pMTHZ/65G47tepLiJSX2paS8lT8v2DQv8OzyDZEWACZpEmY5/APSkbMmYSYiYbmwJk3uho1/d7bcx0vWDJm985k2NqR/BJXYdyH9F4cN6R93X9p87PWJhvt2FIUfmk6ZGZE92pCSoij88MLNTintGjYP49tUSpp1Krf0rbWtwvg2tTvumHK0iqG4o10yAmQEvJYspl0yIhoRiC4ZEY0IJJeM+IniDr5kJG7orEWROyfIPXNZWsqxRv/pZsMJ3We7zvuwaj8DSdr7KlqcKZcPjkG1Czn0jSW6kEP/gGcXcohysCUXctiR6EisNnl2IYf+nTy7kAMshbAkU2PfrfXWd+srWjg7v757fkx6+cVhksbGpOMdv+sTUXVpmKSxMWll+e9vNGpsGBcrpgK+2iIwNrKnYQzq5VRA00YRzBLRLq/QN7vo8gpRKqDo8grRGFRyeYUdicPE6q1nF0Xo68mziyJAdy0labLWoLGlM1ZN1lpEciNu7pRZFJFEKqYiEs4HPNOg6Og3b/eVR7Cvuk355PWyPWwYYf88bcLNR2d1s2GE/X7Aidf6Hj0dJjFATBpY/Ezom3tPhUkMEJMe+33R/NOjy9sws39qR1LE9qwyKjN7E6nlxNsXPrtvuGVTMR/wkaHbvyqeVdLB25ndkd70ySV9QKTXrozIFXumdmWEaK9DdGWEaK9DdGWEvs9QvDJiUvDpbWm2HBEh9umkVyGYmmj6qba9HU07yePO3n+1KH2vWVmHpIkwqf7K0rVL9zbcNqe4Q7GtV6VCKX2jwyRNhElLjz4bM+OPmDBJE9njT6ZJOpxFaxL9eqOWSPTbcZSI2M9VlOjX5PD1/uSZHL6epCiH36n+H/9p3HKA3AvKjiiS+1351jb0dCU+SxoZ072VSk+naHGmDgLOojVJdL0BSSTRlY6cud/390x+XP8RnsmPi3ZNRPLjetv1TH5ctMqLx6A32s1elzl8ptwLLi7bdfdAX0N+saIu0eTRK0dt+dOQ86wox1qs+YNpAx7rYsNu9bsx8z/pUxxd6G63Z5oG8FalKvvD1nUvMGl8njwo5oHepYlESFWypARIIgQmWRKWqiREL1mSZalKxeRTQmCSJeG0FkL0UrqVZzIly1KVhMCk4mSVEL1UnKwSApOKC+GE6KXiNjQhMKkYvAjRSy8HrzgLwYuQqoyDn0PImrDrtThaEAKTLAmLqhGilyzJslSlYoI5ITDJknAaLyF6yZJwYh4hMMmScBovIXrJkmyWqiQEJhWF3gjRS0WpSkJgUnHDgBC9VJSqJAQmFYeEhOill6UqTUcxrEpVDoCfQwhMskc/sFQlIXrJkixLVbIkHFUJgUmWhKMqIXrJkixLVUoP25ga27JUJSEwqXiEkxC9VDzCSQhMKkYwQvTSy1KVpsqDi2SEVCX2J0KejG0s7AX2HcoiBCZZEh71E6KXLMlmqUpCYFJxNEuIXnp5NBtvweIIqcp4+DmEVCW7ho4tzr7De4TAJEvCFkeIXrIkm6UqCYFJRYsjRC+9bHGmhq6BLI6QqsR2QghMKq62EKKX7AQck+yb09k3jrGvLyAEJhV9nBC9ZPfJEqA/uZeqND4vlapkN9sS4YcRSissCS/cEeoCLAkPP4mTxywJqwvYdxqaEJhkSXhST4hesiTLUpUsCZ/yJc7+sCS8PECcC2BJuAsjso9ZEj7lS2RwsiR8ypfIOGNJ+JQvIVXJknBOFrG/zJLwWQVi31B6v7gpyFmWqiQEJhUXpwnRS8XFaUJgUjHLghC9VMyyIAQmFZebCNFLL+enmjoiq1KVuB8jBCZZA8IkQvSSNSBM8qCXNlWHe4FJ4/PkoMEDvUtz6pNtgh6EwCRLwmGQEL1kSTixixCYZEk4NZUQvZTeKG8yAMtSlYTApGJAJUQvFQMqITCpKGBGiF4qCpgRApOKqamE6KViaiohMKmYmkqIXno5NdU0R7AqVYmnGITAJDszxCRC9JJd0sUk+yZQhMAk6x6YRIhesu6BSfYJLhECkywJT1YI0UuWhCUfCIFJloR7H0L0kiXZLFVJCEwq7nMQopeK+xyEwKRipCdELxUjPSEwqShVSYheKkpVEgKTir0PIXrpZckD0xrO3/9iRaoSLwHZt5hkn9AbITDJkvDknxC9ZEk2S1USApOKUpWE6KWXpSpNa31WpSrxUqF9F1vbd+mtBxc9ut8VVbzAz76L2by6zKt4t5F9lzDYJ/Jtn9iuV6UqFSXjPJDdMneBtkk8EVKVipI0HohMmGOBbeIB9h129uAInIlk39EmD46smLoQr0lVmlZsrUpV4u0lQmCSXbXGJEL0ko30mERIVSpuntknkUsITLIr6XijihC9ZEk2S1USApOKWxSE6KXiFgUhMKmY/0GIXno5/8M0IbQqVYnXLezbrvVA3tn0EV6TqhxqIcYRUpV4k9m+7WpCYJJdt8EkQvSSHe9hkn2S24TAJJsIjte3CNFLlmSzVCUhMKl4RogQvVQ8I0QITCqu2xCil15etzEFY6tSlbgvJAQmFU/KEaKX7Hwck+xLSSEEJtmxJSYRopeKYviEwCTb3eL0D0L0kiXZLFVJCEwq7lYSopeKu5WEwKSiVCUheqkoVUkITCpKVRKil16WqjTlaBVDcYeQqsQpXvYlixECk2x3g0mE6CXb3WCSBxdemCrWslQlITCp2PsSopdelqo05fJZlarEqYCEwCSbOY1JhOglu16MSfYlOtp3MQohMMnOjjGJEL1kM9ZGQC+wLFVJCEwqSlUSopeKUpWEwKRiKiAheunlVEDTlpNVqUq8u0YITLKDCUwiRC/ZwQQm2Zc4bN9FQ4TAJKufhEmE6CWrn4RJhMAkuz2JZ6yE6CVLwmvThMAkS8L5gIToJUuyWaqSEJhUHGETopeKI2xCYFJRwIsQvVSUqiQEJhVn9oTopWI+ICEwqZgPSIheejkf0FR5VqUqsY8TApPs8jcmEaKX7PI3JhECk4oHOwjRS3ZXG5M8uFLO1ESWpSoJgUnFHQpC9FJxh4IQmFSUqiREL70sVWma7sNZNCFViVeA7DtKZN/Vg4TAJNtJYBIhesl2EmnQCyxLVRICk4o9HSF66WWpSlNXY1WqEo+IPDhy5n7fX3pFpYlECEyyex2YRIhesvI5mEQITLJrs5hEiF6ym6r4kJ9lqUpCYFJRl4gQvVSUqiQEJhV3qwnRS+ludd3/9an2SWnplJ86n+qR6qP977//wBfubqVwuKvwmWHPrUvLHuHgCzd3FQ6Y5Oga80gvULiFq/CLK0Nvz2uSBgq3dBXeX6ps5KUlo0HhVq7ClZ/e3qjH4+mgcGtX4at+gys8vSkVFG7jKvxhtTt1Lm/oDgq3dRV+v/DW2ilvoHeOcBW+1fbs7FGjkkHhdq7Cg0/+MqTo4fwwvnCkq/C19V8UhJVuBsjtXYVvjx+ePfzhIaBwlKvwyrMJn7YvPBAUjnYVjj2QEjvnShIoHOMqnDUhNPxS1ZdB4Q6uwmk+fhUvR40ChTu6Ck+ufzfiYI/hoHAnV+EijWZ9kRM8EhSOdRXOnX4983z/FFC4s6vw1N6/7WhdD9lGf1fh5LgeUS2ikQ/GuQoX5C4PeiujJyg8wFV4zNWlE589jIx/oKtw/NoZEavqoHeOdxUufGdDqYhNyLsTXIW3vuZ7JaIHeo1EV+Gnx9ab8GJD5IODXIXXN07+8dci6J2TXIWbJuw999U5ZPyDXYV/rhs2rmkC8sEhrsLdMndsmj/uBUBOdhVu+H1B4pXXB4PCQ12F/Rukt9mfNQAUHqar55yaZ/IHIR8c7ircduqCKe8sQW41wlW4TvVT0RuPI09J0VndU6OuRWwaBgqnugpPuL1maYMdKJinuQp/sXZYwIViyAfTXYVbjDz/y5Q1Lzv+B81nuaA=";

{m0, m1, m2} = 
  CoefficientArrays[Last /@ ode, 
   Through[{rZ, rY, rA, rB, rC, rD, rE, rF, rG, rH, rI, rJ, rK, rL, rM, rN, rO, rP, rQ, 
      rR, rS, rT, ra, rb, rc, rd, re, rf, rg, rh, ri, rj, rk, rl, rm, rn, ro, rp, rq, 
      rr, rs, rt}[t$111993]]];
icvalue = Last /@ ic;

Clear@rhs; rhs[var_?VectorQ] := m0 + m1.var + m2.var.var

sol = NDSolveValue[{var'[t] == rhs@var[t], var[0] == icvalue}, 
    var, {t, 0, 10^4}]; // AbsoluteTiming
(* {0.917296, Null} *)    
sol // ListLinePlot

Mathematica graphics

Improve this answer
$\endgroup$
2
  • $\begingroup$ Thanks. Almost any approach works when the number of parameters is small. It is the things like that: github.com/kkkmail/ClmFSharp/blob/master/Clm/Model/Models/… that I am after. And before it can be run, it must be generated. Current F# algorithm (please, see below) needs about 7-9 GB of ram and about 3-5 minutes to generate the file like that, though it can be improved to take much less memory and time. $\endgroup$
    – Konstantin Konstantinov
    Dec 29, 2018 at 0:11
  • $\begingroup$ @KonstantinKonstantinov Well, I know nothing about F#, but since SparseArray is a built-in function, I believe similar performance can be achieved with Mathematica. If you can show us the model with Mathematica code or traditional math notation, I can have a try. $\endgroup$
    – xzczd
    Dec 29, 2018 at 3:24
2
$\begingroup$

I have not been able to determine why it does not work, But there are other methods that return results. One of them is Method -> {"DAEInitialization" -> {"Collocation", "CollocationDirection" -> "Forward"}}

Improve this answer
$\endgroup$
1
$\begingroup$

Technically, this is not an answer, but it works. I moved these models to F# sparse matrices based generator / engine along with a very thin wrapper around ALGLIB (http://www.alglib.net/) vector ODE solver and it just works.

F# takes care of simple, structured, and controlled code generation and ALGLIB is robust enough to handle at least up to 1M variables (yes I tested that and it survived).

Improve this answer
$\endgroup$
3
  • $\begingroup$ If you plug-in the initial conditions found/used with this approach do you get better results with Mathematica? $\endgroup$
    – Anton Antonov
    Dec 28, 2018 at 15:13
  • $\begingroup$ The results are comparable when both work. However, F# / Alglib works for large models where Mathematica fails. My guess is that the main reason is that generation stage gets very demanding even in F# - it could consume 50-60 GB of ram and may run for over an hour. This is due to very large sparse matrices which have only 0.01 to 0.1% of non-zero elements. The resulting F# code file is only about 2 to 10 MB in size and it needs about 100-300 MB to run. $\endgroup$
    – Konstantin Konstantinov
    Dec 28, 2018 at 17:50
  • $\begingroup$ Very interesting! Please put that (or similar comments) in your answer. (+1) $\endgroup$
    – Anton Antonov
    Dec 28, 2018 at 22:28

Your Answer

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

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