numeric-ode (empty) → 0.0.0.0
raw patch · 22 files changed
+1086/−0 lines, 22 filesdep +Chartdep +Chart-cairodep +JuicyPixelssetup-changed
Dependencies added: Chart, Chart-cairo, JuicyPixels, ad, base, colour, data-accessor, data-default-class, diagrams-cairo, diagrams-lib, diagrams-rasterific, foldl, lens, linear, mtl, mwc-probability, mwc-random, numeric-ode, numhask, parallel, parsec, plots, primitive, protolude, reflection, tdigest, template-haskell, text, vector, vector-space
Files
- CHANGES.md +1/−0
- LICENSE +30/−0
- README.md +6/−0
- Setup.hs +2/−0
- diagrams/src_Math_Integrators_StormerVerlet_jupiterOrbit.svg +1/−0
- diagrams/src_Math_Integrators_StormerVerlet_mySquare.svg +1/−0
- numeric-ode.cabal +133/−0
- src/Examples/KeplerProblem.hs +115/−0
- src/Math/Integrators.hs +35/−0
- src/Math/Integrators/ExplicitEuler.hs +18/−0
- src/Math/Integrators/Implicit.hs +64/−0
- src/Math/Integrators/ImplicitEuler.hs +16/−0
- src/Math/Integrators/ImplicitMidpointRule.hs +19/−0
- src/Math/Integrators/Internal.hs +10/−0
- src/Math/Integrators/RK.hs +92/−0
- src/Math/Integrators/RK/Internal.hs +14/−0
- src/Math/Integrators/RK/Parser.hs +71/−0
- src/Math/Integrators/RK/Template.hs +169/−0
- src/Math/Integrators/RK/Types.hs +6/−0
- src/Math/Integrators/StormerVerlet.hs +203/−0
- src/Math/Integrators/StormerVerletAlt.hs +47/−0
- src/Math/Integrators/SympleticEuler.hs +33/−0
+ CHANGES.md view
@@ -0,0 +1,1 @@+TBD
+ LICENSE view
@@ -0,0 +1,30 @@+Copyright (c) 2012, Alexander V Vershilov++All rights reserved.++Redistribution and use in source and binary forms, with or without+modification, are permitted provided that the following conditions are met:++ * Redistributions of source code must retain the above copyright+ notice, this list of conditions and the following disclaimer.++ * Redistributions in binary form must reproduce the above+ copyright notice, this list of conditions and the following+ disclaimer in the documentation and/or other materials provided+ with the distribution.++ * Neither the name of Alexander V Vershilov nor the names of other+ contributors may be used to endorse or promote products derived+ from this software without specific prior written permission.++THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS+"AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT+LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR+A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT+OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,+SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT+LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,+DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY+THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT+(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE+OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+ README.md view
@@ -0,0 +1,6 @@+Small project for different ODE solvers for haskell, in particular+symplectic solvers.++This is very experimental and will change. The Störmer-Verlet+generates a correct orbit for Jupiter but no guarantees are given for+any of the other methods.
+ Setup.hs view
@@ -0,0 +1,2 @@+import Distribution.Simple+main = defaultMain
+ diagrams/src_Math_Integrators_StormerVerlet_jupiterOrbit.svg view
@@ -0,0 +1,1 @@+<svg xmlns="http://www.w3.org/2000/svg" height="400.0000" stroke-opacity="1" viewBox="0 0 500 400" font-size="1" width="500.0000" xmlns:xlink="http://www.w3.org/1999/xlink" stroke="rgb(0,0,0)" version="1.1"><defs></defs><g stroke-linejoin="miter" stroke-opacity="1.0" fill-opacity="1.0" stroke="rgb(0,0,0)" stroke-width="0.0" fill="rgb(255,255,240)" stroke-linecap="butt" stroke-miterlimit="10.0"><path d="M 444.1406,400.0000 l -0.0000,-400.0000 h -388.2812 l -0.0000,400.0000 Z"/></g><g stroke-linejoin="miter" stroke-opacity="1.0" fill-opacity="0.0" stroke="rgb(0,0,0)" stroke-width="0.8" fill="rgb(0,0,0)" stroke-linecap="butt" stroke-miterlimit="10.0"><path d="M 96.4844,352.2383 h 312.5000 M 96.4844,277.0616 h 312.5000 M 96.4844,201.8848 h 312.5000 M 96.4844,126.7080 h 312.5000 M 96.4844,51.5313 h 312.5000 "/></g><g stroke-linejoin="miter" stroke-opacity="1.0" fill-opacity="0.0" stroke="rgb(0,0,0)" stroke-width="1.5763882770434448" fill="rgb(0,0,0)" stroke-linecap="square" stroke-miterlimit="10.0"><path d="M 96.4844,359.3750 v -312.5000 M 408.9844,359.3750 v -312.5000 "/></g><g stroke-linejoin="miter" stroke-opacity="1.0" fill-opacity="0.0" stroke="rgb(0,0,0)" stroke-width="0.4" fill="rgb(0,0,0)" stroke-linecap="butt" stroke-miterlimit="10.0"><path d="M 408.9844,339.7089 h 2.3438 M 408.9844,327.1794 h 2.3438 M 408.9844,314.6499 h 2.3438 M 408.9844,302.1205 h 2.3438 M 408.9844,289.5910 h 2.3438 M 408.9844,264.5321 h 2.3438 M 408.9844,252.0026 h 2.3438 M 408.9844,239.4732 h 2.3438 M 408.9844,226.9437 h 2.3438 M 408.9844,214.4143 h 2.3438 M 408.9844,189.3553 h 2.3438 M 408.9844,176.8259 h 2.3438 M 408.9844,164.2964 h 2.3438 M 408.9844,151.7670 h 2.3438 M 408.9844,139.2375 h 2.3438 M 408.9844,114.1786 h 2.3438 M 408.9844,101.6491 h 2.3438 M 408.9844,89.1197 h 2.3438 M 408.9844,76.5902 h 2.3438 M 408.9844,64.0608 h 2.3438 "/></g><g stroke-linejoin="miter" stroke-opacity="1.0" fill-opacity="0.0" stroke="rgb(0,0,0)" stroke-width="0.4" fill="rgb(0,0,0)" stroke-linecap="butt" stroke-miterlimit="10.0"><path d="M 408.9844,352.2383 h 3.9062 M 408.9844,277.0616 h 3.9062 M 408.9844,201.8848 h 3.9062 M 408.9844,126.7080 h 3.9062 M 408.9844,51.5313 h 3.9062 "/></g><g stroke-linejoin="miter" stroke-opacity="1.0" fill-opacity="0.0" stroke="rgb(0,0,0)" stroke-width="0.4" fill="rgb(0,0,0)" stroke-linecap="butt" stroke-miterlimit="10.0"><path d="M 94.1406,339.7089 h 2.3438 M 94.1406,327.1794 h 2.3438 M 94.1406,314.6499 h 2.3438 M 94.1406,302.1205 h 2.3438 M 94.1406,289.5910 h 2.3438 M 94.1406,264.5321 h 2.3438 M 94.1406,252.0026 h 2.3438 M 94.1406,239.4732 h 2.3438 M 94.1406,226.9437 h 2.3438 M 94.1406,214.4143 h 2.3438 M 94.1406,189.3553 h 2.3438 M 94.1406,176.8259 h 2.3438 M 94.1406,164.2964 h 2.3438 M 94.1406,151.7670 h 2.3438 M 94.1406,139.2375 h 2.3438 M 94.1406,114.1786 h 2.3438 M 94.1406,101.6491 h 2.3438 M 94.1406,89.1197 h 2.3438 M 94.1406,76.5902 h 2.3438 M 94.1406,64.0608 h 2.3438 "/></g><g stroke-linejoin="miter" stroke-opacity="1.0" fill-opacity="0.0" stroke="rgb(0,0,0)" stroke-width="0.4" fill="rgb(0,0,0)" stroke-linecap="butt" stroke-miterlimit="10.0"><path d="M 92.5781,352.2383 h 3.9062 M 92.5781,277.0616 h 3.9062 M 92.5781,201.8848 h 3.9062 M 92.5781,126.7080 h 3.9062 M 92.5781,51.5313 h 3.9062 "/></g><g stroke-linejoin="miter" stroke-opacity="1.0" fill-opacity="1.0" font-size="11.0px" stroke="rgb(0,0,0)" stroke-width="1.5763882770434448" fill="rgb(0,0,0)" stroke-linecap="butt" stroke-miterlimit="10.0"><text dominant-baseline="middle" transform="matrix(1.0000,0.0000,0.0000,1.0000,87.1094,51.5313)" stroke="none" text-anchor="end">1.0e12</text></g><g stroke-linejoin="miter" stroke-opacity="1.0" fill-opacity="1.0" font-size="11.0px" stroke="rgb(0,0,0)" stroke-width="1.5763882770434448" fill="rgb(0,0,0)" stroke-linecap="butt" stroke-miterlimit="10.0"><text dominant-baseline="middle" transform="matrix(1.0000,0.0000,0.0000,1.0000,87.1094,126.7080)" stroke="none" text-anchor="end">5.0e11</text></g><g stroke-linejoin="miter" stroke-opacity="1.0" fill-opacity="1.0" font-size="11.0px" stroke="rgb(0,0,0)" stroke-width="1.5763882770434448" fill="rgb(0,0,0)" stroke-linecap="butt" stroke-miterlimit="10.0"><text dominant-baseline="middle" transform="matrix(1.0000,0.0000,0.0000,1.0000,87.1094,201.8848)" stroke="none" text-anchor="end">0.0</text></g><g stroke-linejoin="miter" stroke-opacity="1.0" fill-opacity="1.0" font-size="11.0px" stroke="rgb(0,0,0)" stroke-width="1.5763882770434448" fill="rgb(0,0,0)" stroke-linecap="butt" stroke-miterlimit="10.0"><text dominant-baseline="middle" transform="matrix(1.0000,0.0000,0.0000,1.0000,87.1094,277.0616)" stroke="none" text-anchor="end">-5.0e11</text></g><g stroke-linejoin="miter" stroke-opacity="1.0" fill-opacity="1.0" font-size="11.0px" stroke="rgb(0,0,0)" stroke-width="1.5763882770434448" fill="rgb(0,0,0)" stroke-linecap="butt" stroke-miterlimit="10.0"><text dominant-baseline="middle" transform="matrix(1.0000,0.0000,0.0000,1.0000,87.1094,352.2383)" stroke="none" text-anchor="end">-1.0e12</text></g><g stroke-linejoin="miter" stroke-opacity="1.0" fill-opacity="0.0" stroke="rgb(0,0,0)" stroke-width="0.8" fill="rgb(0,0,0)" stroke-linecap="butt" stroke-miterlimit="10.0"><path d="M 98.5889,359.3750 v -312.5000 M 173.6617,359.3750 v -312.5000 M 248.7345,359.3750 v -312.5000 M 323.8073,359.3750 v -312.5000 M 398.8802,359.3750 v -312.5000 "/></g><g stroke-linejoin="miter" stroke-opacity="1.0" fill-opacity="0.0" stroke="rgb(0,0,0)" stroke-width="1.5763882770434448" fill="rgb(0,0,0)" stroke-linecap="square" stroke-miterlimit="10.0"><path d="M 96.4844,359.3750 h 312.5000 M 96.4844,46.8750 h 312.5000 "/></g><g stroke-linejoin="miter" stroke-opacity="1.0" fill-opacity="0.0" stroke="rgb(0,0,0)" stroke-width="0.4" fill="rgb(0,0,0)" stroke-linecap="butt" stroke-miterlimit="10.0"><path d="M 111.1010,46.8750 v -2.3438 M 123.6131,46.8750 v -2.3438 M 136.1253,46.8750 v -2.3438 M 148.6374,46.8750 v -2.3438 M 161.1495,46.8750 v -2.3438 M 186.1738,46.8750 v -2.3438 M 198.6860,46.8750 v -2.3438 M 211.1981,46.8750 v -2.3438 M 223.7102,46.8750 v -2.3438 M 236.2224,46.8750 v -2.3438 M 261.2466,46.8750 v -2.3438 M 273.7588,46.8750 v -2.3438 M 286.2709,46.8750 v -2.3438 M 298.7831,46.8750 v -2.3438 M 311.2952,46.8750 v -2.3438 M 336.3195,46.8750 v -2.3438 M 348.8316,46.8750 v -2.3438 M 361.3437,46.8750 v -2.3438 M 373.8559,46.8750 v -2.3438 M 386.3680,46.8750 v -2.3438 "/></g><g stroke-linejoin="miter" stroke-opacity="1.0" fill-opacity="0.0" stroke="rgb(0,0,0)" stroke-width="0.4" fill="rgb(0,0,0)" stroke-linecap="butt" stroke-miterlimit="10.0"><path d="M 98.5889,46.8750 v -3.9062 M 173.6617,46.8750 v -3.9062 M 248.7345,46.8750 v -3.9062 M 323.8073,46.8750 v -3.9062 M 398.8802,46.8750 v -3.9062 "/></g><g stroke-linejoin="miter" stroke-opacity="1.0" fill-opacity="0.0" stroke="rgb(0,0,0)" stroke-width="0.4" fill="rgb(0,0,0)" stroke-linecap="butt" stroke-miterlimit="10.0"><path d="M 111.1010,361.7188 v -2.3438 M 123.6131,361.7188 v -2.3438 M 136.1253,361.7188 v -2.3438 M 148.6374,361.7188 v -2.3438 M 161.1495,361.7188 v -2.3438 M 186.1738,361.7188 v -2.3438 M 198.6860,361.7188 v -2.3438 M 211.1981,361.7188 v -2.3438 M 223.7102,361.7188 v -2.3438 M 236.2224,361.7188 v -2.3438 M 261.2466,361.7188 v -2.3438 M 273.7588,361.7188 v -2.3438 M 286.2709,361.7188 v -2.3438 M 298.7831,361.7188 v -2.3438 M 311.2952,361.7188 v -2.3438 M 336.3195,361.7188 v -2.3438 M 348.8316,361.7188 v -2.3438 M 361.3437,361.7188 v -2.3438 M 373.8559,361.7188 v -2.3438 M 386.3680,361.7188 v -2.3438 "/></g><g stroke-linejoin="miter" stroke-opacity="1.0" fill-opacity="0.0" stroke="rgb(0,0,0)" stroke-width="0.4" fill="rgb(0,0,0)" stroke-linecap="butt" stroke-miterlimit="10.0"><path d="M 98.5889,363.2812 v -3.9062 M 173.6617,363.2812 v -3.9062 M 248.7345,363.2812 v -3.9062 M 323.8073,363.2812 v -3.9062 M 398.8802,363.2812 v -3.9062 "/></g><g stroke-linejoin="miter" stroke-opacity="1.0" fill-opacity="1.0" font-size="11.0px" stroke="rgb(0,0,0)" stroke-width="1.5763882770434448" fill="rgb(0,0,0)" stroke-linecap="butt" stroke-miterlimit="10.0"><text dominant-baseline="text-before-edge" transform="matrix(1.0000,0.0000,0.0000,1.0000,398.8802,368.7500)" stroke="none" text-anchor="middle">1.0e12</text></g><g stroke-linejoin="miter" stroke-opacity="1.0" fill-opacity="1.0" font-size="11.0px" stroke="rgb(0,0,0)" stroke-width="1.5763882770434448" fill="rgb(0,0,0)" stroke-linecap="butt" stroke-miterlimit="10.0"><text dominant-baseline="text-before-edge" transform="matrix(1.0000,0.0000,0.0000,1.0000,323.8073,368.7500)" stroke="none" text-anchor="middle">5.0e11</text></g><g stroke-linejoin="miter" stroke-opacity="1.0" fill-opacity="1.0" font-size="11.0px" stroke="rgb(0,0,0)" stroke-width="1.5763882770434448" fill="rgb(0,0,0)" stroke-linecap="butt" stroke-miterlimit="10.0"><text dominant-baseline="text-before-edge" transform="matrix(1.0000,0.0000,0.0000,1.0000,248.7345,368.7500)" stroke="none" text-anchor="middle">0.0</text></g><g stroke-linejoin="miter" stroke-opacity="1.0" fill-opacity="1.0" font-size="11.0px" stroke="rgb(0,0,0)" stroke-width="1.5763882770434448" fill="rgb(0,0,0)" stroke-linecap="butt" stroke-miterlimit="10.0"><text dominant-baseline="text-before-edge" transform="matrix(1.0000,0.0000,0.0000,1.0000,173.6617,368.7500)" stroke="none" text-anchor="middle">-5.0e11</text></g><g stroke-linejoin="miter" stroke-opacity="1.0" fill-opacity="1.0" font-size="11.0px" stroke="rgb(0,0,0)" stroke-width="1.5763882770434448" fill="rgb(0,0,0)" stroke-linecap="butt" stroke-miterlimit="10.0"><text dominant-baseline="text-before-edge" transform="matrix(1.0000,0.0000,0.0000,1.0000,98.5889,368.7500)" stroke="none" text-anchor="middle">-1.0e12</text></g><g stroke-linejoin="miter" stroke-opacity="1.0" fill-opacity="0.0" stroke="rgb(255,127,0)" stroke-width="3.0" fill="rgb(0,0,0)" stroke-linecap="butt" stroke-miterlimit="10.0"><defs><clipPath id="myClip1"><path d="M 408.9844,359.3750 l -0.0000,-312.5000 h -312.5000 l -0.0000,312.5000 Z"/></clipPath></defs><g clip-path="url(#myClip1)"><path d="M 137.5406,201.8848 l 1.2141,17.8100 l 3.8817,17.3780 l 6.4307,16.5326 l 8.7942,15.3060 l 10.9148,13.7417 l 12.7460,11.8914 l 14.2535,9.8121 l 15.4147,7.5633 l 16.2185,5.2037 l 16.6635,2.7902 l 16.7568,0.3757 l 16.5125,-1.9915 l 15.9499,-4.2684 l 15.0924,-6.4170 l 13.9664,-8.4048 l 12.6000,-10.2038 l 11.0229,-11.7905 l 9.2652,-13.1454 l 7.3579,-14.2524 l 5.3319,-15.0987 l 3.2187,-15.6743 l 1.0500,-15.9717 l -1.1421,-15.9861 l -3.3250,-15.7152 l -5.4654,-15.1595 l -7.5292,-14.3223 l -9.4818,-13.2101 l -11.2877,-11.8334 l -12.9111,-10.2066 l -14.3162,-8.3494 l -15.4676,-6.2867 l -16.3317,-4.0495 l -16.8774,-1.6756 l -17.0778,0.7907 l -16.9118,3.2985 l -16.3660,5.7909 l -15.4366,8.2061 l -14.1313,10.4795 l -12.4702,12.5461 l -10.4867,14.3432 l -8.2273,15.8141 l -5.7503,16.9113 l -3.1234,17.5995 "/></g></g><g stroke-linejoin="miter" stroke-opacity="1.0" fill-opacity="1.0" opacity="0.55" stroke="rgb(0,0,0)" stroke-width="1.5763882770434448" fill="rgb(165,42,42)" stroke-linecap="butt" stroke-miterlimit="10.0"><defs><clipPath id="myClip2"><path d="M 408.9844,359.3750 l -0.0000,-312.5000 h -312.5000 l -0.0000,312.5000 Z"/></clipPath></defs><g clip-path="url(#myClip2)"><path d="M 170.1958,141.7313 c 0.0000,-8.3038 -6.7222,-15.0354 -15.0146 -15.0354c -8.2923,-0.0000 -15.0146,6.7316 -15.0146 15.0354c -0.0000,8.3038 6.7222,15.0354 15.0146 15.0354c 8.2923,0.0000 15.0146,-6.7316 15.0146 -15.0354Z"/></g></g><g stroke-linejoin="miter" stroke-opacity="1.0" fill-opacity="1.0" opacity="0.5" stroke="rgb(0,0,0)" stroke-width="1.5763882770434448" fill="rgb(165,42,42)" stroke-linecap="butt" stroke-miterlimit="10.0"><defs><clipPath id="myClip3"><path d="M 408.9844,359.3750 l -0.0000,-312.5000 h -312.5000 l -0.0000,312.5000 Z"/></clipPath></defs><g clip-path="url(#myClip3)"><path d="M 222.7206,96.1564 c 0.0000,-8.3038 -6.7222,-15.0354 -15.0146 -15.0354c -8.2923,-0.0000 -15.0146,6.7316 -15.0146 15.0354c -0.0000,8.3038 6.7222,15.0354 15.0146 15.0354c 8.2923,0.0000 15.0146,-6.7316 15.0146 -15.0354Z"/></g></g><g stroke-linejoin="miter" stroke-opacity="1.0" fill-opacity="1.0" opacity="0.45" stroke="rgb(0,0,0)" stroke-width="1.5763882770434448" fill="rgb(165,42,42)" stroke-linecap="butt" stroke-miterlimit="10.0"><defs><clipPath id="myClip4"><path d="M 408.9844,359.3750 l -0.0000,-312.5000 h -312.5000 l -0.0000,312.5000 Z"/></clipPath></defs><g clip-path="url(#myClip4)"><path d="M 289.9535,87.9520 c 0.0000,-8.3038 -6.7222,-15.0354 -15.0146 -15.0354c -8.2923,-0.0000 -15.0146,6.7316 -15.0146 15.0354c -0.0000,8.3038 6.7222,15.0354 15.0146 15.0354c 8.2923,0.0000 15.0146,-6.7316 15.0146 -15.0354Z"/></g></g><g stroke-linejoin="miter" stroke-opacity="1.0" fill-opacity="1.0" opacity="0.4" stroke="rgb(0,0,0)" stroke-width="1.5763882770434448" fill="rgb(165,42,42)" stroke-linecap="butt" stroke-miterlimit="10.0"><defs><clipPath id="myClip5"><path d="M 408.9844,359.3750 l -0.0000,-312.5000 h -312.5000 l -0.0000,312.5000 Z"/></clipPath></defs><g clip-path="url(#myClip5)"><path d="M 348.9802,116.8442 c 0.0000,-8.3038 -6.7222,-15.0354 -15.0146 -15.0354c -8.2923,-0.0000 -15.0146,6.7316 -15.0146 15.0354c -0.0000,8.3038 6.7222,15.0354 15.0146 15.0354c 8.2923,0.0000 15.0146,-6.7316 15.0146 -15.0354Z"/></g></g><g stroke-linejoin="miter" stroke-opacity="1.0" fill-opacity="1.0" opacity="0.35000000000000003" stroke="rgb(0,0,0)" stroke-width="1.5763882770434448" fill="rgb(165,42,42)" stroke-linecap="butt" stroke-miterlimit="10.0"><defs><clipPath id="myClip6"><path d="M 408.9844,359.3750 l -0.0000,-312.5000 h -312.5000 l -0.0000,312.5000 Z"/></clipPath></defs><g clip-path="url(#myClip6)"><path d="M 382.7443,171.3696 c 0.0000,-8.3038 -6.7222,-15.0354 -15.0146 -15.0354c -8.2923,-0.0000 -15.0146,6.7316 -15.0146 15.0354c -0.0000,8.3038 6.7222,15.0354 15.0146 15.0354c 8.2923,0.0000 15.0146,-6.7316 15.0146 -15.0354Z"/></g></g><g stroke-linejoin="miter" stroke-opacity="1.0" fill-opacity="1.0" opacity="0.30000000000000004" stroke="rgb(0,0,0)" stroke-width="1.5763882770434448" fill="rgb(165,42,42)" stroke-linecap="butt" stroke-miterlimit="10.0"><defs><clipPath id="myClip7"><path d="M 408.9844,359.3750 l -0.0000,-312.5000 h -312.5000 l -0.0000,312.5000 Z"/></clipPath></defs><g clip-path="url(#myClip7)"><path d="M 382.9427,234.7168 c 0.0000,-8.3038 -6.7222,-15.0354 -15.0146 -15.0354c -8.2923,-0.0000 -15.0146,6.7316 -15.0146 15.0354c -0.0000,8.3038 6.7222,15.0354 15.0146 15.0354c 8.2923,0.0000 15.0146,-6.7316 15.0146 -15.0354Z"/></g></g><g stroke-linejoin="miter" stroke-opacity="1.0" fill-opacity="1.0" opacity="0.25000000000000006" stroke="rgb(0,0,0)" stroke-width="1.5763882770434448" fill="rgb(165,42,42)" stroke-linecap="butt" stroke-miterlimit="10.0"><defs><clipPath id="myClip8"><path d="M 408.9844,359.3750 l -0.0000,-312.5000 h -312.5000 l -0.0000,312.5000 Z"/></clipPath></defs><g clip-path="url(#myClip8)"><path d="M 349.9648,289.0039 c 0.0000,-8.3038 -6.7222,-15.0354 -15.0146 -15.0354c -8.2923,-0.0000 -15.0146,6.7316 -15.0146 15.0354c -0.0000,8.3038 6.7222,15.0354 15.0146 15.0354c 8.2923,0.0000 15.0146,-6.7316 15.0146 -15.0354Z"/></g></g><g stroke-linejoin="miter" stroke-opacity="1.0" fill-opacity="1.0" opacity="0.20000000000000004" stroke="rgb(0,0,0)" stroke-width="1.5763882770434448" fill="rgb(165,42,42)" stroke-linecap="butt" stroke-miterlimit="10.0"><defs><clipPath id="myClip9"><path d="M 408.9844,359.3750 l -0.0000,-312.5000 h -312.5000 l -0.0000,312.5000 Z"/></clipPath></defs><g clip-path="url(#myClip9)"><path d="M 292.3561,318.2980 c 0.0000,-8.3038 -6.7222,-15.0354 -15.0146 -15.0354c -8.2923,-0.0000 -15.0146,6.7316 -15.0146 15.0354c -0.0000,8.3038 6.7222,15.0354 15.0146 15.0354c 8.2923,0.0000 15.0146,-6.7316 15.0146 -15.0354Z"/></g></g><g stroke-linejoin="miter" stroke-opacity="1.0" fill-opacity="1.0" opacity="0.15000000000000002" stroke="rgb(0,0,0)" stroke-width="1.5763882770434448" fill="rgb(165,42,42)" stroke-linecap="butt" stroke-miterlimit="10.0"><defs><clipPath id="myClip10"><path d="M 408.9844,359.3750 l -0.0000,-312.5000 h -312.5000 l -0.0000,312.5000 Z"/></clipPath></defs><g clip-path="url(#myClip10)"><path d="M 226.2048,311.9199 c 0.0000,-8.3038 -6.7222,-15.0354 -15.0146 -15.0354c -8.2923,-0.0000 -15.0146,6.7316 -15.0146 15.0354c -0.0000,8.3038 6.7222,15.0354 15.0146 15.0354c 8.2923,0.0000 15.0146,-6.7316 15.0146 -15.0354Z"/></g></g><g stroke-linejoin="miter" stroke-opacity="1.0" fill-opacity="1.0" opacity="0.1" stroke="rgb(0,0,0)" stroke-width="1.5763882770434448" fill="rgb(165,42,42)" stroke-linecap="butt" stroke-miterlimit="10.0"><defs><clipPath id="myClip11"><path d="M 408.9844,359.3750 l -0.0000,-312.5000 h -312.5000 l -0.0000,312.5000 Z"/></clipPath></defs><g clip-path="url(#myClip11)"><path d="M 172.8758,268.9115 c 0.0000,-8.3038 -6.7222,-15.0354 -15.0146 -15.0354c -8.2923,-0.0000 -15.0146,6.7316 -15.0146 15.0354c -0.0000,8.3038 6.7222,15.0354 15.0146 15.0354c 8.2923,0.0000 15.0146,-6.7316 15.0146 -15.0354Z"/></g></g><g stroke-linejoin="miter" stroke-opacity="1.0" fill-opacity="1.0" opacity="5.0e-2" stroke="rgb(0,0,0)" stroke-width="1.5763882770434448" fill="rgb(165,42,42)" stroke-linecap="butt" stroke-miterlimit="10.0"><defs><clipPath id="myClip12"><path d="M 408.9844,359.3750 l -0.0000,-312.5000 h -312.5000 l -0.0000,312.5000 Z"/></clipPath></defs><g clip-path="url(#myClip12)"><path d="M 152.5552,201.8848 c 0.0000,-8.3038 -6.7222,-15.0354 -15.0146 -15.0354c -8.2923,-0.0000 -15.0146,6.7316 -15.0146 15.0354c -0.0000,8.3038 6.7222,15.0354 15.0146 15.0354c 8.2923,0.0000 15.0146,-6.7316 15.0146 -15.0354Z"/></g></g><g stroke-linejoin="miter" stroke-opacity="1.0" fill-opacity="1.0" stroke="rgb(0,0,0)" stroke-width="1.5763882770434448" fill="rgb(255,255,0)" stroke-linecap="butt" stroke-miterlimit="10.0"><defs><clipPath id="myClip13"><path d="M 408.9844,359.3750 l -0.0000,-312.5000 h -312.5000 l -0.0000,312.5000 Z"/></clipPath></defs><g clip-path="url(#myClip13)"><path d="M 263.7491,201.8848 c 0.0000,-8.3038 -6.7222,-15.0354 -15.0146 -15.0354c -8.2923,-0.0000 -15.0146,6.7316 -15.0146 15.0354c -0.0000,8.3038 6.7222,15.0354 15.0146 15.0354c 8.2923,0.0000 15.0146,-6.7316 15.0146 -15.0354Z"/></g></g><g stroke-linejoin="miter" stroke-opacity="1.0" fill-opacity="1.0" font-size="11.0px" stroke="rgb(0,0,0)" stroke-width="1.5763882770434448" fill="rgb(0,0,0)" stroke-linecap="butt" stroke-miterlimit="10.0"><text dominant-baseline="text-after-edge" transform="matrix(1.0000,0.0000,0.0000,1.0000,252.7344,31.2500)" stroke="none" text-anchor="middle"></text></g></svg>
+ diagrams/src_Math_Integrators_StormerVerlet_mySquare.svg view
@@ -0,0 +1,1 @@+<svg xmlns="http://www.w3.org/2000/svg" height="300.0000" stroke-opacity="1" viewBox="0 0 200 300" font-size="1" width="200.0000" xmlns:xlink="http://www.w3.org/1999/xlink" stroke="rgb(0,0,0)" version="1.1"><defs></defs><g fill-opacity="1.0" fill="rgb(0,0,255)"><path d="M 200.0000,211.6809 l -38.3191,-161.6809 l -161.6809,38.3191 l 38.3191,161.6809 Z"/></g></svg>
+ numeric-ode.cabal view
@@ -0,0 +1,133 @@+name: numeric-ode+version: 0.0.0.0+synopsis: Ode solvers+description: Some ode solvers, e.g., Störmer-Verlet+homepage: https://github.com/qnikst/numeric-ode+license: BSD3+license-file: LICENSE+author: Alexander V Vershilov, Dominic Steinitz+maintainer: dominic@steinitz.org+copyright: Alexander V Vershilov, Dominic Steinitz+category: Math+build-type: Simple+cabal-version: >=1.18+extra-source-files: README.md, CHANGES.md, diagrams/*.svg+extra-doc-files: diagrams/*.svg++source-repository head+ type: git+ location: https://github.com/qnikst/numeric-ode++library+ default-language: Haskell2010+ hs-source-dirs: src+ exposed-modules:+ Math.Integrators+ Math.Integrators.ExplicitEuler+ Math.Integrators.ImplicitEuler+ Math.Integrators.ImplicitMidpointRule+ Math.Integrators.SympleticEuler+ Math.Integrators.StormerVerlet+ Math.Integrators.StormerVerletAlt+ Math.Integrators.RK+ Math.Integrators.Implicit+ Math.Integrators.Internal+ Math.Integrators.RK.Internal+ Math.Integrators.RK.Parser+ Math.Integrators.RK.Template+ Math.Integrators.RK.Types++ ghc-options: -Wall+ -- other-modules:+ build-depends: base>=4 && < 5,+ vector>=0.9 && <1.1,+ parallel>=3.2 && <3.3,+ parsec == 3.1.*,+ template-haskell,+ linear,+ lens,+ primitive>=0.4 && <0.7,+ text,+ protolude,+ mwc-random,+ mwc-probability,+ primitive,+ ad,+ reflection,+ tdigest,+ -- chart-unit,+ numhask,+ foldl++ other-extensions: TypeFamilies+ FlexibleContexts+ BangPatterns+ QuasiQuotes++executable Kepler+ hs-source-dirs: src/Examples+ main-is: KeplerProblem.hs+ ghc-options:+ build-depends: base,+ numeric-ode,+ vector>=0.9 && <1.0,+ vector-space>=0.8 && <0.11,+ colour,+ linear,+ data-default-class,+ diagrams-lib,+ diagrams-cairo,+ Chart,+ Chart-cairo,+ data-accessor,+ diagrams-rasterific,+ diagrams-lib,+ JuicyPixels,+ plots,+ mtl++ default-language: Haskell2010++-- executable TestChart+-- hs-source-dirs: src/Examples+-- main-is: TestChart.hs+-- ghc-options:+-- build-depends: base >= 4.7 && < 5,+-- chart-unit,+-- protolude,+-- foldl,+-- text,+-- numhask,+-- -- for data examples+-- mwc-random,+-- mwc-probability,+-- primitive,+-- ad,+-- reflection,+-- tdigest,+-- diagrams-cairo,+-- diagrams-lib,+-- JuicyPixels+-- default-language: Haskell2010++-- executable TestRasterific+-- hs-source-dirs: src/Examples+-- main-is: TestRasterific.hs+-- ghc-options:+-- build-depends: base >= 4.7 && < 5,+-- chart-unit,+-- protolude,+-- foldl,+-- text,+-- numhask,+-- -- for data examples+-- mwc-random,+-- mwc-probability,+-- primitive,+-- ad,+-- reflection,+-- tdigest,+-- diagrams-rasterific,+-- diagrams-lib,+-- JuicyPixels+-- default-language: Haskell2010
+ src/Examples/KeplerProblem.hs view
@@ -0,0 +1,115 @@+{-# LANGUAGE NegativeLiterals #-}+{-# LANGUAGE TypeFamilies #-}+{-# LANGUAGE FlexibleContexts #-}++{-# OPTIONS_GHC -Wall #-}++module Main (main) where++import qualified Data.Vector as V+import Control.Monad.ST++import Math.Integrators.StormerVerlet+import Math.Integrators++import qualified Linear as L+import Linear.V+import Data.Maybe ( fromJust )++import Diagrams.Prelude+import Diagrams.Backend.CmdLine+import Diagrams.Backend.Rasterific.CmdLine++import Control.Monad+import Control.Monad.State.Class++import Plots+++gConst :: Double+gConst = 6.67384e-11++nStepsTwoPlanets :: Int+nStepsTwoPlanets = 44++stepTwoPlanets :: Double+stepTwoPlanets = 24 * 60 * 60 * 100++sunMass, jupiterMass :: Double+sunMass = 1.9889e30+jupiterMass = 1.8986e27++jupiterPerihelion :: Double+jupiterPerihelion = 7.405736e11++jupiterV :: [Double]+jupiterV = [-1.0965244901087316e02, -1.3710001990210707e04, 0.0]++jupiterQ :: [Double]+jupiterQ = [negate jupiterPerihelion, 0.0, 0.0]++sunV :: [Double]+sunV = [0.0, 0.0, 0.0]++sunQ :: [Double]+sunQ = [0.0, 0.0, 0.0]++tm :: V.Vector Double+tm = V.enumFromStepN 0 stepTwoPlanets nStepsTwoPlanets++kepler :: L.V2 (L.V3 Double) -> L.V2 (L.V3 Double)+kepler (L.V2 q1 q2) =+ let r = q2 L.^-^ q1+ ri = r `L.dot` r+ rr = ri * (sqrt ri)+ q1' = pure gConst * r / pure rr+ q2' = negate q1'+ q1'' = q1' * pure sunMass+ q2'' = q2' * pure jupiterMass+ in L.V2 q1'' q2''++listToV3 :: [a] -> L.V3 a+listToV3 [x, y, z] = fromV . fromJust . fromVector . V.fromList $ [x, y, z]+listToV3 xs = error $ "Only supply 3 elements not: " ++ show (length xs)++initPQs :: L.V2 (L.V2 (L.V3 Double))+initPQs = L.V2 (L.V2 (listToV3 jupiterV) (listToV3 sunV))+ (L.V2 (listToV3 jupiterQ) (listToV3 sunQ))++result1 :: V.Vector (L.V2 (L.V2 (L.V3 Double)))+result1 = runST $ integrateV (\h -> stormerVerlet2 kepler (pure h)) initPQs tm++preMorePts :: [(Double, Double)]+preMorePts = map (\(L.V2 _ (L.V2 (L.V3 x y _z) _)) -> (x,y)) (V.toList result1)++morePts :: [P2 Double]+morePts = map p2 $ preMorePts++addPoint :: (Plotable (Diagram B) b, MonadState (Axis b V2 Double) m) =>+ Double -> (Double, Double) -> m ()+addPoint o (x, y) = addPlotable'+ ((circle 1e11 :: Diagram B) #+ fc brown #+ opacity o #+ translate (r2 (x, y)))++jSaxis :: Axis B V2 Double+jSaxis = r2Axis &~ do+ addPlotable' ((circle 1e11 :: Diagram B) # fc yellow)+ let l = length preMorePts+ let os = [0.05,0.1..]+ let ps = take (l `div` 4) [0,4..]+ zipWithM_ addPoint os (map (preMorePts!!) ps)+ linePlot' $ map unp2 $ take 200 morePts++displayHeader :: FilePath -> Diagram B -> IO ()+displayHeader fn =+ mainRender ( DiagramOpts (Just 900) (Just 700) fn+ , DiagramLoopOpts False Nothing 0+ )++main :: IO ()+main = do+ displayHeader "other/jupiter-sun-line.png" (renderAxis jSaxis # bg ivory)+ putStrLn "Finished"+
+ src/Math/Integrators.hs view
@@ -0,0 +1,35 @@+-- | Math integrators if a high level module for different ODE integrators+-- This module provides high-level wrappers over different integration methods+--+module Math.Integrators + where++import Data.Vector (Vector,(!))+import Data.Vector.Mutable+import Control.Monad.Primitive+import qualified Data.Vector as V++import Math.Integrators.Internal++{-|+ Integrate ODE equation using fixed steps set by a vector, and returns a vector+ of solutions corrensdonded to times that was requested.+ It takes Vector of time points as a parameter and returns a vector of results+ -}+integrateV :: PrimMonad m => Integrator a -- ^ Internal integrator+ -> a -- ^ initial value+ -> Vector Double -- ^ vector of time points+ -> m (Vector a) -- ^ vector of solution+integrateV integrator initial times = do+ out <- new (V.length times) + write out 0 initial+ compute initial 1 out+ V.unsafeFreeze out+ where+ compute y i out | i == V.length times = return () + | otherwise = do+ let h = (times ! i) - (times ! (i-1))+ y' = integrator h y+ write out i y'+ compute y' (i+1) out+
+ src/Math/Integrators/ExplicitEuler.hs view
@@ -0,0 +1,18 @@+{-# LANGUAGE FlexibleContexts #-}+-- |+-- Module: Math.Integrators.ExplicitEuler.+--+-- Basic integrator using Euler method. It has following properies:+-- * allows to solve systems of the first order+-- * this method is not symplectic and tends to loose energy+--+module Math.Integrators.ExplicitEuler+ where++import Linear++-- | Integrator of form+--+-- \[ \Phi[h] : y_{n+1} = y_n + h f (y_n) \]+explicitEuler :: (Num (f a), Floating a, Additive f) => (f a -> f a) -> a -> f a -> f a+explicitEuler f = \h y -> y ^+^ h *^ (f y)
+ src/Math/Integrators/Implicit.hs view
@@ -0,0 +1,64 @@+{-# LANGUAGE FlexibleContexts #-}+-- | Helpers for implicit integration methods+--+-- TODO: add possibility to make function to create initial value+-- TODO: add possibility to break on step+-- TODO: add possibility to add different initial value based+-- on y0, f+-- TODO: add seq-pseq to make this stuff strict+-- TODO: add Newton iterations+module Math.Integrators.Implicit+ ( -- * types+ ImplicitSolver+ -- * solvers+ , fixedPointSolver+ , fixedPoint+ -- * helpers+ , breakNormR+ , breakNormIR+ )+ where++import Linear+import Control.Lens++-- | Implicit solver type+type ImplicitSolver a = (a -> a) -- ^ implicit method+ -> (Int -> a -> a -> Bool) -- ^ breakRule+ -> a -- ^ initial value+ -> a -- ^ final value++-- | Fixed point method it iterates function f until it will break "" will+-- be reached then it returns one but last iteration+--+fixedPointSolver :: ImplicitSolver a+fixedPointSolver f break' y0 = inner 0 y0+ where + inner i y = let y' = f y+ i' = i+1+ in if break' i y y'+ then y'+ else inner i' y'++fixedPoint :: (a -> a) -- ^ function+ -> (a -> a -> Bool) -- ^ break rule+ -> a -- ^ initial value+ -> a -- ^ result+fixedPoint f break' y0 = + let y1 = f y0+ in if break' y0 y1+ then y0+ else fixedPoint f break' y1++-- | simple break rule that will break evaluatioin when value less then Eps+breakNormR :: Double -> Double -> Bool+breakNormR eps y = abs y < eps++-- | same as @breakNormR@ but assume that inner type is an +-- instance of InnerField, so it's possible to use innerproduct to find norm+-- N.B function uses $||v||^2 < eps$, so epsilon should be pre evaluated+breakNormIR :: (Metric f, Floating a, Ord a, Num (f a)) => f a -> a -> Bool+breakNormIR v eps = quadrance v < eps+++
+ src/Math/Integrators/ImplicitEuler.hs view
@@ -0,0 +1,16 @@+module Math.Integrators.ImplicitEuler+ ( implicitEuler+ ) where++import Linear++import Math.Integrators.Implicit+import Math.Integrators.Internal++eps :: Floating a => a+eps = 1e-14++implicitEuler :: (Metric f, Ord a, Additive f, Num (f a), Floating a)+ => (f a -> f a) -> a -> f a -> f a+implicitEuler f = \h y ->+ fixedPoint (\x -> y ^+^ (h *^ (f x))) (\x1 x2 -> breakNormIR (x1^-^x2) eps) y
+ src/Math/Integrators/ImplicitMidpointRule.hs view
@@ -0,0 +1,19 @@+{-# LANGUAGE FlexibleContexts #-}+module Math.Integrators.ImplicitMidpointRule+ ( imr+ ) where++import Linear++import Math.Integrators.Implicit+import Math.Integrators.Internal++eps :: Floating a => a+eps = 1e-14++imr :: (Metric f, Num (f a), Floating a, Ord a)+ => (f a -> f a) -> a -> f a -> f a+imr f = \h y ->+ fixedPoint (\x -> y ^+^ h *^ ( f ( (y^+^x)^/2) ))+ (\x1 x2 -> breakNormIR (x1 ^-^ x2) eps)+ y
+ src/Math/Integrators/Internal.hs view
@@ -0,0 +1,10 @@+module Math.Integrators.Internal+ where++{- | Integrator function+ - \Phi [h] |-> y_0 -> y_1+ -}+type Integrator a = Double -- ^ Step+ -> a -- ^ Initial value+ -> a -- ^ Next value+
+ src/Math/Integrators/RK.hs view
@@ -0,0 +1,92 @@+{-# LANGUAGE QuasiQuotes, FlexibleContexts #-}+{-# OPTIONS_GHC -Wwarn #-}+-- | Runge-Kutta module +-- TODO: add description and history notes+-- add informations about methods properties+module Math.Integrators.RK+ ( -- * explicit methods+-- rk45+-- , rk46+-- -- * implicit methods+-- , gauss4+-- , gauss6+-- , lobattoIIIA4+-- , lobattoIIIA6+-- , lobattoIIIB4+ ) where++-- import Linear++-- import Math.Integrators.RK.Template+-- import Math.Integrators.RK.Types+-- import Math.Integrators.Internal+-- import Math.Integrators.Implicit+++{-+rk45 :: (VectorSpace a, Floating (Scalar a)) => (Double -> a -> a) -> Integrator (Double,a)+rk45 = [qrk|+0 |+0.5 | 0.5+0.5 | 0 & 0.5+1 | 0 & 0 & 1+- - + - - - - - - - - + | 1/6 & 2/6 & 2/6 & 1/6+|]++rk46 :: (VectorSpace a, Floating (Scalar a)) => (Double -> a -> a) -> Integrator (Double,a)+rk46 = [qrk|+0 |+1/3 | 1/3+2/3 | -1/3 & 1+1 | 1 & -1 & 1+- - + - - - - - - - - + | 1/8 & 3/8 & 3/8 & 1/8+|]+++gauss4 :: (VectorSpace a, Floating (Scalar a)) => (ImplicitRkType (a,a)) -> (Double -> a -> a) -> Integrator (Double,a)+gauss4 = [qrk|+0.5 - sqrt(3)/6 | 0.25 & 0.25 - sqrt(3)/6+0.5 + sqrt(3)/6 | 0.25 + sqrt(3)/6 & 1/4+- - - - - - - - + - - - - - - - - - - - + | 0.5 & 0.5+|]++gauss6 :: (VectorSpace a, Floating (Scalar a)) => (ImplicitRkType (a,a,a)) -> (Double -> a -> a) -> Integrator (Double,a)+gauss6 = [qrk|+0.5 - sqrt(15)/10 | 5/36 & 2/9 - sqrt(15)/15 & 5/36 - sqrt(15)/30+0.5 | 5/36 + sqrt(15)/24 & 2/9 & 5/36 - sqrt(15)/24+0.5 + sqrt(15)/10 | 5/36 + sqrt(15)/30 & 2/9+sqrt(15)/15 & 5/36+- - - - - - - - - + - - - - - - - - - - - - - - - - - - - - - - -+ | 5/18 & 4/9 & 5/18+|]++lobattoIIIA4 :: (VectorSpace a, Floating (Scalar a)) => (ImplicitRkType (a,a,a)) -> (Double -> a -> a) -> Integrator (Double,a)+lobattoIIIA4 = [qrk|+0 | 0 & 0 & 0+0.5 | 5/24 & 1/3 & -1/24+1 | 1/6 & 2/3 & 1/6+- - + - - - - - - - - - - + | 1/6 & 2/3 & 1/6+|]++lobattoIIIA6 :: (VectorSpace a, Floating (Scalar a)) => (ImplicitRkType (a,a,a,a)) -> (Double -> a -> a) -> Integrator (Double,a)+lobattoIIIA6 = [qrk|+0 | 0 & 0 & 0 & 0+(5-sqrt(5))/10 | (11+sqrt(5))/120 & (25-sqrt(5))/120 & (25 - 13 *sqrt(5)/120) & (-1+sqrt(5))/120+(5+sqrt(5))/10 | (11-sqrt(5))/120 & (25+13*sqrt(5))/120 & (25+sqrt(5))/120 & (-1-sqrt(5))/120+1 | 1/12 & 5/12 & 5/12 & 1/12+- - - - - - - - + - - - -+ | 1/12 & 5/12 & 5/12 & 1/12+|]++lobattoIIIB4 :: (VectorSpace a, Floating (Scalar a)) => (ImplicitRkType (a,a,a)) -> (Double -> a -> a) -> Integrator (Double,a)+lobattoIIIB4 = [qrk|+0 | 1/6 & -1/6 & 0+0.5 | 1/6 & 1/3 & 0+1 | 1/6 & 5/6 & 0+- - + - - - - - - - - -+ | 1/6 & 2/3 & 1/6+|]+-}
+ src/Math/Integrators/RK/Internal.hs view
@@ -0,0 +1,14 @@+module Math.Integrators.RK.Internal+ ( MExp(..)+ , isExplicit+ )+ where+++-- | Internal type that users by+data MExp = Delimeter | Row (Maybe Double,[Double]) deriving (Show,Eq)++isExplicit :: [MExp] -> Bool+isExplicit = (all (\(i,(Row (_,x))) -> i> length x)) . (zip [1..]) . top+ where + top = takeWhile (/= Delimeter)
+ src/Math/Integrators/RK/Parser.hs view
@@ -0,0 +1,71 @@+{-# OPTIONS_GHC -Wwarn #-} -- We need this option, because we want to remove this module in future+module Math.Integrators.RK.Parser+ ( readMatrixTable+ )+ where++import Data.Maybe++-- Parsec stuff+import Text.Parsec+import qualified Text.Parsec.Token as P+import Text.Parsec.Language (haskellDef) +import Text.Parsec.Expr+import Text.Parsec.String++import Math.Integrators.RK.Internal++readMatrixTable :: String -> [MExp]+readMatrixTable = + mapMaybe go . lines+ where+ go ('-':s) = Just Delimeter+ go ('#':s) = Nothing+ go (ls) | null.filter (==' ')$ ls = Nothing+ | otherwise = + let (lhs,_:rhs) = span (/='|') ls+ l = case filter (/= ' ') lhs of+ "" -> Nothing+ ls -> Just $! erun expr ls+ in Just $ Row (l, map (erun expr) $! grp '&' rhs)+ grp c s = case dropWhile (==c) s of+ "" -> []+ s' -> if any (/=' ') w then w : grp c c'' else grp c c''+ where (w,c'') = break (==c) s'++lexer = P.makeTokenParser haskellDef { P.reservedOpNames = ["*","/","+","-","sqrt","sin","cos"] }++whiteSpace= P.whiteSpace lexer+lexeme = P.lexeme lexer+symbol = P.symbol lexer+float = P.float lexer+parens = P.parens lexer+natural = P.natural lexer+identifier= P.identifier lexer+reserved = P.reserved lexer+reservedOp= P.reservedOp lexer++expr :: Parser Double+expr = buildExpressionParser table factor+ <?> "expression"+factor = parens expr+ <|> try float+ <|> fmap realToFrac natural+ <?> "simple expression"+table = [ [prefix "-" negate]+ , [prefix "sqrt" sqrt,prefix "sin" sin,prefix "cos" cos]+ , [op "*" (*) AssocLeft, op "/" (/) AssocLeft]+ , [op "+" (+) AssocLeft, op "-" (-) AssocLeft]+ ] + where+ op s f assoc = Infix (do{ reservedOp s; return f} <?> "operator") assoc+ prefix s f = Prefix (do { reservedOp s; return f} <?> "prefix")+++erun :: Parser Double -> String -> Double+erun p input = erun' (do { whiteSpace ; x <- p ; eof; return x})+ where+ erun' p' = case (parse p' "" input) of+ Left err -> error $ "Parse error at "++(show err)+ Right x -> x+
+ src/Math/Integrators/RK/Template.hs view
@@ -0,0 +1,169 @@+{-# LANGUAGE TemplateHaskell #-}+{-# OPTIONS_GHC -Wwarn #-} -- this module will be removed in future versions+module Math.Integrators.RK.Template+ where++import Data.Maybe++import Language.Haskell.TH.Quote+import Language.Haskell.TH++import Control.Monad++import Math.Integrators.RK.Types+import Math.Integrators.RK.Internal+import Math.Integrators.RK.Parser+++qrk :: QuasiQuoter+qrk = QuasiQuoter {quoteExp = x}+ where+ x s = rk $! readMatrixTable s++-----------------------------------------------------------+-- List of helpers+jv = Just . VarE+jv' = Just . varE+jld = Just . LitE. {-DoublePrimL-} RationalL . toRational+ld = LitE . RationalL . toRational+plus = VarE $! mkName "+"+vplus = VarE $! mkName "^+^"+vmult = VarE $! mkName "*^"+mult = VarE $! mkName "*"+ld' = litE . RationalL . toRational+plus' = varE $! mkName "+"+vplus'= varE $! mkName "^+^"+vmult'= varE $! mkName "*^"+mult' = varE $! mkName "*"+vN s = varE (mkName s)++foldOp op = foldl1 (\x y -> infixE (Just x) op (Just y))++realToFracN = varE (mkName "realToFrac")+zeroVN = varE (mkName "zeroV")+f = mkName "f"+t = mkName "t"+h = mkName "h"+y = mkName "y"+tpy = mkName "tpy"++rk :: [MExp] -> Q Exp+rk mExp = do+ let (ab,_:c:[]) = break (==Delimeter) mExp+ lenA = length ab+ kn <- forM [1..lenA] (\_ -> newName "k")+ let kvv = zip kn ab+ ks' <- forM kvv $ \(k,r) -> do+ t <- rkt1 r kn+ return $ ValD (VarP k) (NormalB t) []+ y' <- rkt2 c kn+ if isExplicit mExp+ then return $ LamE [VarP f, VarP h, TupP [VarP t,VarP y]] $ LetE ks' y'+ else irk mExp+ where + rkt1 (Row (Just c,ls)) ks = do+ let ft = infixE (jv' t) plus' (Just $ infixE (Just $ ld' c) mult' (jv' h))+ st = if null ls + then (varE y)+ else+ infixE (jv' y) vplus'+ (Just $ infixE (Just $ appE realToFracN (varE h)) vmult'+ (Just $ foldOp vplus' $+ zipWith (\k l -> infixE (Just $ appE realToFracN (ld' l)) vmult' (jv' k)) + ks + ls+ )+ )+ appE (appE (varE f) ft) st+ rkt2 (Row (_,ls)) ks = do+ tupE [ infixE (jv' t) plus' (jv' h)+ , infixE (jv' y) vplus'+ (Just $ infixE (Just $ appE realToFracN (varE h)) vmult'+ (Just $ foldOp vplus' $+ zipWith (\k l -> infixE (Just $ appE realToFracN (ld' l)) vmult' (jv' k)) ks ls+ )+ )+ ]+++test = [Row (Just 1,[2,3]),Row (Just 4,[5,6]),Delimeter, Row (Nothing, [7,8])]++irk mExp = do+ fpoint' <- fpoint mExp+ lamE [varP tpy, varP f, varP h, tupP [varP t,varP y]] $ + caseE (varE tpy) [match (conP (mkName "Math.Integrators.RK.Types.FixedPoint") [varP $ mkName "breakRule"])+ (normalB (fpointRun mExp)) + fpoint'+ ,match (conP (mkName "Math.Integrators.RK.Types.NewtonIteration") []) (normalB (varE $ mkName "undefined")) []+ ]++fpointRun mExp = do+ let (ab,_:(Row (_,ls)):[]) = break (==Delimeter) mExp+ lenA = length ab+ zs <- forM [1..lenA] (\_ -> newName "z")+ letE [valD (tupP $ map varP zs)+ (normalB $ + appE + (appE + (appE + (varE $! mkName "Math.Integrators.Implicit.fixedPointSolver") + (varE $! mkName "method")+ ) + (varE $ mkName "breakRule")+ ) + (tupE $ replicate lenA zeroVN) {- TODO: give avaliability to user -}+ ) + []]+ (appE (varE (mkName "solution")) (tupE $ map varE zs))++fpoint mExp = do+ let (ab,_:(Row (_,ls)):[]) = break (==Delimeter) mExp+ lenA = length ab+ zs <- forM [1..lenA] (\_ -> newName "z")+ zs' <- forM [1..lenA] (const $ newName "z'")++ return $ + [ funD (mkName "method") [clause [tupP $ map varP zs] (normalB $ letE (map (topRow zs) $! zip zs' ab) (tupE $ map varE zs')) [] ]+ , funD (mkName "solution") [clause [tupP $ map varP zs] (normalB $ solutionRow zs ls) []]+ ]+ where+ topRow zs (x,(Row (Just c,ls))) = + valD (varP x) + (normalB $ infixE (Just $ appE realToFracN (varE h)) + vmult'+ (Just $ foldOp vplus' $+ zipWith (\z l -> infixE (Just $ appE realToFracN (ld' l))+ vmult'+ (Just $ appE (appE (varE f) + (infixE (jv' t) + plus' + (Just $ infixE (Just $ appE realToFracN $ varE h) mult' (Just $ ld' c))+ )+ )+ (infixE (jv' y) vplus' (jv' z))+ )+ )+ zs+ ls+ )+ )+ []+ topRow _ _ = error "not a row"+ solutionRow zs ls = + (tupE [infixE (jv' t) plus' (jv' h)+ ,infixE (jv' y) + vplus' + (Just $ infixE (Just $ appE realToFracN $ varE h) + vmult'+ (Just $ foldOp vplus' + (zipWith (\z b -> infixE (Just $ appE realToFracN + (ld' b))+ vmult'+ (jv' z))+ zs+ ls+ ) + )+ )+ ]+ )
+ src/Math/Integrators/RK/Types.hs view
@@ -0,0 +1,6 @@+module Math.Integrators.RK.Types+ where+++-- | type implicit solver+data ImplicitRkType a = FixedPoint (Int -> a -> a -> Bool) | NewtonIteration
+ src/Math/Integrators/StormerVerlet.hs view
@@ -0,0 +1,203 @@+{-# LANGUAGE FlexibleContexts #-}++-- |+-- Module: Math.Integrators.StormerVerlet+--+--+-- Störmer-Verlet is an order 2 symplectic method. This means it will+-- preserve the Hamiltonian for the system the differential equations+-- describe, for example, important for modelling planetary motion;+-- the application of something like the much-loved Runge-Kutta 4th+-- order method would either model the planet spiralling toward or+-- away from the Sun!+--+-- Here's a diagram showing the orbit of Jupiter around the Sun.+-- +-- <<diagrams/src_Math_Integrators_StormerVerlet_jupiterOrbit.svg#diagram=jupiterOrbit&height=400&width=500>>+--+-- To create this, consider the \(n\)-body problem. The Hamiltonian is+--+-- \[+-- {\mathbb H} = \frac{1}{2}\sum_{i=0}^n \frac{p_i^\top p_i}{m_i} - \frac{G}{2}\sum_{i=0}^n\sum_{j \neq i} \frac{m_i m_j}{\|q_i - q_j\|}+-- \]+--+-- Apply Hamilton's equations will gives \(2n\) first order+-- equations. To use 'stormerVerlet2' this needs to be \(n\) second order+-- equations. In this case, the Lagrangian is easy+--+-- \[+-- {\mathcal{L}} = \frac{1}{2}\sum_{i=0}^n \frac{p_i^\top p_i}{m_i} + \frac{G}{2}\sum_{i=0}^n\sum_{j \neq i} \frac{m_i m_j}{\|q_i - q_j\|}+-- \]+--+-- Applying Lagrange's equation+--+-- \[+-- \frac{\mathrm{d}}{\mathrm{d}t}\bigg(\frac{\partial{\mathcal{L}}}{\partial\dot{q}_j}\bigg) = \frac{\partial{\mathcal{L}}}{\partial{q}_j}+-- \]+--+-- gives+--+-- \[+-- m_j\ddot{q}_j = G\sum_{k \neq j}m_k m_j \frac{q_k - q_j}{\|q_k - q_j\|^3}+-- \]+--+-- For \(n = 2\) this gives+--+-- \[+-- \begin{aligned}+-- \ddot{q}_1 &= m_2G\frac{q_1 - q_2}{\|q_1 - q_2\|^3} \\+-- \ddot{q}_2 &= m_1G\frac{q_2 - q_1}{\|q_2 - q_1\|^3}+-- \end{aligned}+-- \]+--+-- > {-# LANGUAGE NegativeLiterals #-}+-- > {-# LANGUAGE TypeFamilies #-}+-- > {-# LANGUAGE FlexibleContexts #-}+-- > {-# LANGUAGE MultiParamTypeClasses #-}+-- > +-- > import qualified Data.Vector as V+-- > import Control.Monad.ST+-- > +-- > import Math.Integrators.StormerVerlet+-- > import Math.Integrators+-- > +-- > import qualified Linear as L+-- > import Linear.V+-- > import Data.Maybe ( fromJust )+-- > +-- > import Diagrams.Prelude+-- > +-- > import Control.Monad+-- > import Control.Monad.State.Class+-- > +-- > import Plots+-- >+-- > -- First some constants describing the system+-- >+-- > gConst :: Double+-- > gConst = 6.67384e-11+-- > +-- > nStepsTwoPlanets :: Int+-- > nStepsTwoPlanets = 44+-- >+-- > -- A step size of 100 days!+-- >+-- > stepTwoPlanets :: Double+-- > stepTwoPlanets = 24 * 60 * 60 * 100+-- > +-- > sunMass, jupiterMass :: Double+-- > sunMass = 1.9889e30+-- > jupiterMass = 1.8986e27+-- > +-- > jupiterPerihelion :: Double+-- > jupiterPerihelion = 7.405736e11+-- > +-- > jupiterV :: L.V3 Double+-- > jupiterV = L.V3 (-1.0965244901087316e02) (-1.3710001990210707e04) 0.0+-- > +-- > jupiterQ :: L.V3 Double+-- > jupiterQ = L.V3 (-jupiterPerihelion) 0.0 0.0+-- > +-- > sunV :: L.V3 Double+-- > sunV = L.V3 0.0 0.0 0.0+-- > +-- > sunQ :: L.V3 Double+-- > sunQ = L.V3 0.0 0.0 0.0+-- >+-- > -- The right hand side of the second order differential equation system.+-- >+-- > kepler :: L.V2 (L.V3 Double) -> L.V2 (L.V3 Double)+-- > kepler (L.V2 q1 q2) =+-- > let r = q2 L.^-^ q1+-- > ri = r `L.dot` r+-- > rr = ri * (sqrt ri)+-- > q1' = pure gConst * r / pure rr+-- > q2' = negate q1'+-- > q1'' = q1' * pure sunMass+-- > q2'' = q2' * pure jupiterMass+-- > in L.V2 q1'' q2''+-- >+-- > -- Initial values+-- >+-- > initPQs :: L.V2 (L.V2 (L.V3 Double))+-- > initPQs = L.V2 (L.V2 jupiterV sunV) (L.V2 jupiterQ sunQ)+-- >+-- > -- Steps at which to evolve the system+-- >+-- > tm :: V.Vector Double+-- > tm = V.enumFromStepN 0 stepTwoPlanets nStepsTwoPlanets+-- >+-- > -- The results+-- >+-- > result1 :: V.Vector (L.V2 (L.V2 (L.V3 Double)))+-- > result1 = runST $ integrateV (\h -> stormerVerlet2 kepler (pure h)) initPQs tm+-- > +-- > preMorePts :: [(Double, Double)]+-- > preMorePts = map (\(L.V2 _ (L.V2 (L.V3 x y _z) _)) -> (x,y)) (V.toList result1)+-- > +-- > morePts :: [P2 Double]+-- > morePts = map p2 $ preMorePts+-- >+-- > -- Finally plot the results+-- >+-- > addPoint :: (Plotable (Diagram B) b, MonadState (Axis b V2 Double) m) =>+-- > Double -> (Double, Double) -> m ()+-- > addPoint o (x, y) = addPlotable'+-- > ((circle 1e11 :: Diagram B) #+-- > fc brown #+-- > opacity o #+-- > translate (r2 (x, y)))+-- > +-- > jSaxis :: Axis B V2 Double+-- > jSaxis = r2Axis &~ do+-- > addPlotable' ((circle 1e11 :: Diagram B) # fc yellow)+-- > let l = length preMorePts+-- > let os = [0.05,0.1..]+-- > let ps = take (l `div` 4) [0,4..]+-- > zipWithM_ addPoint os (map (preMorePts!!) ps)+-- > linePlot' $ map unp2 $ take 200 morePts+-- > +-- > jupiterOrbit = renderAxis jSaxis # bg ivory+--++module Math.Integrators.StormerVerlet+ ( stormerVerlet2H+ , stormerVerlet2+ ) where++import Linear+import Control.Lens++-- | Störmer-Verlet integration scheme for systems of the form+-- \(\mathbb{H}(p,q) = T(p,q) + V(p,q)\)+stormerVerlet2H :: (Applicative f, Num (f a), Fractional a) =>+ a -- ^ Step size+ -> (f a -> f a) -- ^ \(\frac{\partial H}{\partial q}\)+ -> (f a -> f a) -- ^ \(\frac{\partial H}{\partial p}\)+ -> V2 (f a) -- ^ Current \((p, q)\) as a 2-dimensional vector+ -> V2 (f a) -- ^ New \((p, q)\) as a 2-dimensional vector+stormerVerlet2H hh nablaQ nablaP prev = V2 qNew pNew+ where+ h2 = hh / 2+ hhs = pure hh+ hh2s = pure h2+ qsPrev = prev ^. _x+ psPrev = prev ^. _y+ pp2 = psPrev - hh2s * nablaQ qsPrev+ qNew = qsPrev + hhs * nablaP pp2+ pNew = pp2 - hh2s * nablaQ qNew++-- | Störmer-Verlet integration scheme for system: \(\ddot{\mathbf{q}} = f(\mathbf{q})\)+stormerVerlet2 :: (Applicative f, Num (f a), Fractional a)+ => (f a -> f a) -- ^ \(f\)+ -> a -- ^ Step size+ -> V2 (f a) -- ^ Current \((p, q)\) as a 2-dimensional vector+ -> V2 (f a) -- ^ New \((p, q)\) as a 2-dimensional vector+stormerVerlet2 f h prev =+ let h' = h+ h2' = 0.5 * h+ p1 = prev ^. _x + pure h2' * (f (prev ^. _y))+ q' = prev ^. _y + pure h' * p1+ p' = p1 + pure h2' * (f q')+ in V2 p' q'+
+ src/Math/Integrators/StormerVerletAlt.hs view
@@ -0,0 +1,47 @@+{-# OPTIONS_GHC -Wall #-}+{-# OPTIONS_GHC -fno-warn-type-defaults #-}++module Math.Integrators.StormerVerletAlt where++import Linear+import Control.Lens++oneStepH98 :: (Applicative f, Num (f a), Fractional a) =>+ a -- ^ Step size+ -> (f a -> f a) -- ^ \(\frac{\partial H}{\partial q}\)+ -> (f a -> f a) -- ^ \(\frac{\partial H}{\partial p}\)+ -> V2 (f a) -- ^ Current \((p, q)\) as a 2-dimensional vector+ -> V2 (f a) -- ^ New \((p, q)\) as a 2-dimensional vector+oneStepH98 hh nablaQ nablaP prev = V2 qNew pNew+ where+ h2 = hh / 2+ hhs = pure hh+ hh2s = pure h2+ qsPrev = prev ^. _x+ psPrev = prev ^. _y+ pp2 = psPrev - hh2s * nablaQ qsPrev+ qNew = qsPrev + hhs * nablaP pp2+ pNew = pp2 - hh2s * nablaQ qNew++-- And now can apply this to the two body problem with the following+-- derivatives of the Hamiltonian.++nablaQ' :: V2 Double -> V2 Double+nablaQ' qs = V2 (qq1 / r) (qq2 / r)+ where+ qq1 = qs ^. _x+ qq2 = qs ^. _y+ r = (qq1 ^ 2 + qq2 ^ 2) ** (3/2)++nablaP' :: V2 Double -> V2 Double+nablaP' ps = ps++e, q10, q20, p10, p20 :: Double+e = 0.6+q10 = 1 - e+q20 = 0.0+p10 = 0.0+p20 = sqrt ((1 + e) / (1 - e))++inits :: V2 (V2 Double)+inits = V2 (V2 q10 q20) (V2 p10 p20)
+ src/Math/Integrators/SympleticEuler.hs view
@@ -0,0 +1,33 @@+{-# LANGUAGE FlexibleContexts #-}+module Math.Integrators.SympleticEuler+ where++import Linear+import Control.Lens++import Math.Integrators.Implicit++eps :: Floating a => a+eps = 1e-10++sympleticEuler1 :: (Metric f, Num (f a), Floating a, Ord a)+ => (f a -> f a -> f a)+ -> (f a -> f a -> f a) + -> a -- ^ Step size+ -> V2 (f a) -- ^ Current \((p,q)\) as a 2-dimentional vector+ -> V2 (f a) -- ^ New \((p, q)\) as a 2-dimetional vector+sympleticEuler1 f g = \h prev ->+ -- explicit coordinate+ let u' = (prev^._x) ^+^ h *^ (f (prev^._x) v')+ -- implicit coordinate+ v' = fixedPoint (\x -> (prev^._y) ^+^ h *^ (g (prev^._x) x))+ (\x1 x2 -> breakNormIR (x1^-^x2) eps) (prev^._x)+ in V2 u' v'++{-+sEuler2 :: ((a->a->a),(a->a->a)) -> Double -> (a,a) -> (a,a)+sEuler2 (a,b) h (u,v) =+ let u' = u + ( h * (a u' v) )+ v' = v + ( h * (b u' v) )+ in (u',v')+-}