{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "{-# LANGUAGE RankNTypes, TypeFamilies #-}\n",
    "import Data.VectorSpace\n",
    "import Data.Manifold.Types"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "type ODESolver\n",
    "   = forall v . (VectorSpace v, Scalar v ~ ℝ)\n",
    "       => (v -> v)   -- ^ Function to integrate\n",
    "       -> ℝ          -- ^ Time step\n",
    "       -> v          -- ^ Start state\n",
    "       -> [v]        -- ^ Sequence of calculated steps"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "euler :: ODESolver\n",
    "euler f h y = y : euler f h (y ^+^ h*^f y)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "rk₄ :: ODESolver\n",
    "rk₄ f h y = y : rk₄ f h (y ^+^ (h/6)*^(k₁ ^+^ 2*^k₂ ^+^ 2*^k₃ ^+^ k₄))\n",
    " where k₁ = f y\n",
    "       k₂ = f $ y ^+^ (h/2)*^k₁\n",
    "       k₃ = f $ y ^+^ (h/2)*^k₂\n",
    "       k₄ = f $ y ^+^ h*^k₃"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "GraphWindowSpecR2{lBound=-13.874009118754717, rBound=25.510352576781575, bBound=-19.831199555474484, tBound=35.94014689231086, xResolution=640, yResolution=480}"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import Graphics.Dynamic.Plot.R2\n",
    "import Data.List\n",
    "plotWindow\n",
    "    [ legendName name\n",
    "     $ plotLatest\n",
    "       [ lineSegPlot [(x,y) | (x,y,_) <- segment]\n",
    "       | segment <- take 400 <$> iterate (drop 10) traject ]\n",
    "    | (solver, name) <- [(euler, \"Euler\"), (rk₄,\"RK₄\")]\n",
    "    , let traject = solver (\\(x,y,z) -> (σ*(y-x), x*(ρ-z), x*y - β*z))\n",
    "                           0.002\n",
    "                           (10,11,12)\n",
    "           where ρ=28; σ=10; β=8/3\n",
    "    ]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Haskell",
   "language": "haskell",
   "name": "haskell"
  },
  "language_info": {
   "codemirror_mode": "ihaskell",
   "file_extension": ".hs",
   "name": "haskell",
   "version": "8.2.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}