By reducing tolerances cvsRoberts_dns now agrees with manual for every decade from...
authorilb@NIH.GOV <ilb@NIH.GOV@ba61647d-9d00-f842-95cd-605cb4296b96>
Thu, 8 Mar 2018 22:55:07 +0000 (22:55 +0000)
committerilb@NIH.GOV <ilb@NIH.GOV@ba61647d-9d00-f842-95cd-605cb4296b96>
Thu, 8 Mar 2018 22:55:07 +0000 (22:55 +0000)
git-svn-id: https://citdcbmipav.cit.nih.gov/repos-pub/mipav/trunk@15408 ba61647d-9d00-f842-95cd-605cb4296b96

mipav/src/gov/nih/mipav/model/algorithms/CVODES.java

index 9585e85..e971c19 100644 (file)
@@ -614,20 +614,62 @@ public abstract class CVODES {
                 * tolerance. Output is printed in decades from t = .4 to t = 4.e10.\r
                 * Run statistics (optional outputs) are printed at the end.\r
                 * -----------------------------------------------------------------*/\r
+               //Example manual and MIPAV agree for t = 0.4, 4.0, 400.0, t = 4.0E3, t = 4.0E4, t = 4.0E5,\r
+           // t = 4.0E6, t = 4.0E7, and t = 4.0E8\r
+           // MIPAV stops agreeing with example manual at 4.0E9\r
+               // At t = 4.0E10 the example manual gives:\r
+               // y[0] = 6.5181E-8 y[1] = 2.6072E-13 y[2] = 1.0000\r
+               // but the supplied check_ans routine gives:\r
+               // ref.data[0] = 5.2083495894337328e-08;\r
+               // ref.data[1] = 2.0833399429795671e-13;\r
+               // ref.data[2] = 9.9999994791629776e-01;\r
+               /*Example manual: at t = 0.4 y[0] = 0.98517 y[1] = 3.3864E-5 y[2] = 1.4794E-2\r
+               MIPAV: at t = 0.4 y[0] = 0.9851721138611713 y[1] = 3.3863953789795866E-5 y[2] = 0.014794022184947894\r
+               Example manual: at t = 4.0 y[0] = 0.90552 y[1] = 2.2405E-5 y[2] = 9.4459E-2\r
+               MIPAV: at t = 4.0 y[0] = 0.9055186785921142 y[1] = 2.2404756875816905E-5 y[2] = 0.09445891665808899\r
+               Example manual: at t = 40.0 y[0] = 0.71583 y[1] = 9.1856E-6 y[2] = 0.28416\r
+               MIPAV: at t = 40.0 y[0] = 0.715827068745404 y[1] = 9.185534765344522E-6 y[2] = 0.28416374572813036\r
+               Example manual: at t = 400.0 y[0] = 0.45052 y[1] = 3.2229E-6 y[2] = 0.54947\r
+               MIPAV: at t = 400.0 y[0] = 0.4505186684223784 y[1] = 3.222901440811379E-6 y[2] = 0.5494781087190751\r
+               Example manual: at t = 4.0E3 y[0] = 0.18317 y[1] = 8.9403E-7 y[2] = 0.81683\r
+               MIPAV: at t = 4000.0 y[0] = 0.18320225775817142 y[1] = 8.942371253468589E-7 y[2] = 0.8167968478397172\r
+               Example manual: at t = 4.0E4 y[0] = 3.8977E-2 y[1] = 1.6215E-7 y[2] = 0.96102\r
+               MIPAV: at t = 40000.0 y[0] = 0.038983377128004884 y[1] = 1.6217683145641717E-7 y[2] = 0.9610164625843524\r
+               Example manual: at t = 4.0E5 y[0] = 4.9387E-3 y[1] = 1.9852E-8 y[2] = 0.99506\r
+               MIPAV: at t = 400000.0 y[0] = 0.004938274376153642 y[1] = 1.98499403707653E-8 y[2] = 0.9950617019719638\r
+               Example manual: at t = 4.0E6 y[0] = 5.1684E-4 y[1] = 2.0684E-9 y[2] = 0.99948\r
+               MIPAV: at t = 4000000.0 y[0] = 5.168097052506471E-4 y[1] = 2.0682949105260644E-9 y[2] = 0.9994831816270917\r
+               Example manual: at t = 4.0E7 y[0] = 5.2039E-5 y[1] = 2.0817E-10 y[2] = 0.99995\r
+               MIPAV: at t = 2.0795462081963886E7 y[0] = 1.0E-4 y[1] = 4.0004167668416476E-10 y[2] = 0.9998946217354228\r
+               Roots found are -1 and 0\r
+               Example manual: at t = 4.0E7 y[0] = 5.2039E-5 y[1] = 2.0817E-10 y[2] = 0.99995\r
+               MIPAV: at t = 4.0E7 y[0] = 5.203026008088301E-5 y[1] = 2.0813195051402364E-10 y[2] = 0.9999469745454322\r
+               Example manual at t = 4.0E8 y[0] = 5.2106E-6 y[1] = 2.0842E-11 y[2] = 0.99999\r
+               MIPAV: at t = 4.0E8 y[0] = 5.231316609955509E-6 y[1] = 2.092540715323315E-11 y[2] = 0.9999929133739337\r
+               Example manual at t = 4.0E9 y[0] = 5.1881E-7 y[1] = 2.0752E-12 y[2] = 1.0000\r
+               MIPAV: at t = 4.0E9 y[0] = 7.480696384365109E-7 y[1] = 2.9922793356220927E-12 y[2] = 0.9999980891645214\r
+               Example manual at t = 4.0E10 y[0] = 6.5181E-8 y[1] = 2.6072E-13 y[2] = 1.0000\r
+        MIPAV: at t = 4.0E10 y[0] = 6.181145073087833E-7 y[1] = 2.472445298118015E-12 y[2] = 1.0000048993342436*/\r
+               \r
                \r
                /** Problem Constants */\r
                final int NEQ = 3; // Number of equations\r
                final double Y1 = 1.0; // Initial y components\r
                final double Y2 = 0.0;\r
                final double Y3 = 0.0;\r
-               final double RTOL = 1.0E-4; // scalar relative tolerance\r
-               final double ATOL1 = 1.0E-8; // vector absolute tolerance components\r
-               final double ATOL2 = 1.0E-14;\r
-               final double ATOL3 = 1.0E-6;\r
+               //final double RTOL = 1.0E-4; // scalar relative tolerance\r
+               final double RTOL = 1.0E-12;\r
+               //final double ATOL1 = 1.0E-8; // vector absolute tolerance components\r
+               final double ATOL1 = 1.0E-12;\r
+               //final double ATOL2 = 1.0E-14;\r
+               final double ATOL2 = 1.0E-15;\r
+               //final double ATOL3 = 1.0E-6;\r
+               final double ATOL3 = 1.0E-12;\r
                final double T0 = 0.0; // initial time\r
                final double T1 = 0.4; // first output time\r
                final double TMULT = 10.0; // output time factor\r
-               final int NOUT = 12; // number of output times\r
+               //final int NOUT = 12; // number of output times\r
+               final int NOUT = 12;\r
                int f = cvsRoberts_dns;\r
                int g = cvsRoberts_dns;\r
                int Jac = cvsRoberts_dns;\r
@@ -730,7 +772,45 @@ public abstract class CVODES {
                \r
                while (true) {\r
                        flag = CVode(cvode_mem, tout, y, t, CV_NORMAL);\r
-                       System.out.println("At t = " + t[0] + " y[0] = " + y.data[0] + " y[1] = " + y.data[1] + " y[2] = " + y.data[2]);\r
+                       switch (iout) {\r
+                       case 0:\r
+                               System.out.println("Example manual: at t = 0.4 y[0] = 0.98517 y[1] = 3.3864E-5 y[2] = 1.4794E-2");\r
+                               break;\r
+                       case 1:\r
+                               System.out.println("Example manual: at t = 4.0 y[0] = 0.90552 y[1] = 2.2405E-5 y[2] = 9.4459E-2");\r
+                               break;\r
+                       case 2:\r
+                               System.out.println("Example manual: at t = 40.0 y[0] = 0.71583 y[1] = 9.1856E-6 y[2] = 0.28416");\r
+                               break;\r
+                       case 3:\r
+                               System.out.println("Example manual: at t = 400.0 y[0] = 0.45052 y[1] = 3.2229E-6 y[2] = 0.54947");\r
+                               break;\r
+                       case 4:\r
+                               System.out.println("Example manual: at t = 4.0E3 y[0] = 0.18317 y[1] = 8.9403E-7 y[2] = 0.81683");\r
+                               break;\r
+                       case 5:\r
+                               System.out.println("Example manual: at t = 4.0E4 y[0] = 3.8977E-2 y[1] = 1.6215E-7 y[2] = 0.96102");\r
+                               break;\r
+                       case 6:\r
+                               System.out.println("Example manual: at t = 4.0E5 y[0] = 4.9387E-3 y[1] = 1.9852E-8 y[2] = 0.99506");\r
+                               break;\r
+                       case 7:\r
+                               System.out.println("Example manual: at t = 4.0E6 y[0] = 5.1684E-4 y[1] = 2.0684E-9 y[2] = 0.99948");\r
+                               break;\r
+                       case 8:\r
+                               System.out.println("Example manual: at t = 4.0E7 y[0] = 5.2039E-5 y[1] = 2.0817E-10 y[2] = 0.99995");\r
+                               break;\r
+                       case 9:\r
+                               System.out.println("Example manual at t = 4.0E8 y[0] = 5.2106E-6 y[1] = 2.0842E-11 y[2] = 0.99999");\r
+                               break;\r
+                       case 10:\r
+                               System.out.println("Example manual at t = 4.0E9 y[0] = 5.1881E-7 y[1] = 2.0752E-12 y[2] = 1.0000");\r
+                               break;\r
+                       case 11:\r
+                               System.out.println("Example manual at t = 4.0E10 y[0] = 6.5181E-8 y[1] = 2.6072E-13 y[2] = 1.0000");\r
+                               break;\r
+                       }\r
+                       System.out.println("MIPAV: at t = " + t[0] + " y[0] = " + y.data[0] + " y[1] = " + y.data[1] + " y[2] = " + y.data[2]);\r
                        \r
                        if (flag == CV_ROOT_RETURN) {\r
                            flagr = CVodeGetRootInfo(cvode_mem, rootsfound)     ;\r
@@ -794,7 +874,7 @@ public abstract class CVODES {
        }\r
        \r
        /* compare the solution at the final time 4e10s to a reference solution computed\r
-          using a relative tolerance of 1e-8 and absoltue tolerance of 1e-14 */\r
+          using a relative tolerance of 1e-8 and absolute tolerance of 1e-14 */\r
        private int check_ans(NVector y, double rtol, NVector atol)\r
        {\r
          int      passfail=0;        /* answer pass (0) or fail (1) flag */  \r