J. Cash BVP problems

j_cash_bench.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Some problems from J. Cash suite"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import numpy as np"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "from scipy.integrate import solve_bvp"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 43,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "from scipy.special import erf, gamma, airy"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 44,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "pi = np.pi"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 45,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def check_accuracy(y, y_true):\n",
    "    d = (y - y_true) / np.maximum(1, np.abs(y_true))\n",
    "    return np.mean(np.abs(d))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 46,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "class TestProblem(object):\n",
    "    def __init__(self, a, b, va, vb, eps):\n",
    "        self.x = np.linspace(a, b, 5)\n",
    "        self.y = np.zeros((2, 5))\n",
    "        self.va = va\n",
    "        self.vb = vb\n",
    "        self.eps = eps\n",
    "        \n",
    "    def rhs(self, x, y):\n",
    "        raise NotImplementedError\n",
    "        \n",
    "    def ref(self, x):\n",
    "        raise NotImplementedError\n",
    "        \n",
    "    def fun(self, x, y):\n",
    "        return np.vstack((\n",
    "            y[1],\n",
    "            self.rhs(x, y) / self.eps\n",
    "        ))\n",
    "        \n",
    "    def bc(self, ya, yb):\n",
    "        return np.array([ya[0] - self.va, yb[0] - self.vb])\n",
    "    \n",
    "    def run(self, tol=1e-3):\n",
    "        res = solve_bvp(self.fun, self.bc, self.x, self.y, tol=tol)\n",
    "        return res.success, res.x.size, check_accuracy(res.y[0], self.ref(res.x))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 47,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "class Problem4(TestProblem):\n",
    "    def __init__(self, eps):\n",
    "        super(Problem4, self).__init__(-1, 1, 1 + np.exp(-2), 1 + np.exp(-2 * (1 + 1 / eps)), eps)\n",
    "                                        \n",
    "    def rhs(self, x, y):\n",
    "        return -y[1] + (1 + self.eps) * y[0]\n",
    "    \n",
    "    def ref(self, x):\n",
    "        return np.exp(x - 1) + np.exp(-(1 + eps) * (1 + x) / eps)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 48,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "class Problem6(TestProblem): \n",
    "    def __init__(self, eps):\n",
    "        super(Problem6, self).__init__(-1, 1, -2, 0, eps)\n",
    "    \n",
    "    def rhs(self, x, y):\n",
    "        return -(x * y[1] + self.eps * pi**2 * np.cos(pi * x) + pi * x * np.sin(pi * x))\n",
    "            \n",
    "    def ref(self, x):\n",
    "        k = (2 * self.eps) ** 0.5\n",
    "        return np.cos(pi * x) + erf(x / k) / erf(1 / k) "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 49,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "class Problem7(TestProblem):\n",
    "    def __init__(self, eps):\n",
    "        super(Problem7, self).__init__(-1, 1, -1, 1, eps)\n",
    "        \n",
    "    def rhs(self, x, y):\n",
    "        eps = self.eps\n",
    "        return -x * y[1] + y[0] - (1 + eps * pi**2) * np.cos(pi * x) - pi * x * np.sin(pi * x)\n",
    "    \n",
    "    def ref(self, x):\n",
    "        k = (2 * self.eps) ** 0.5\n",
    "        eps = self.eps\n",
    "        return (np.cos(pi * x) + x + \n",
    "                (x * erf(x / k) + k / pi**0.5 * np.exp(-0.5 * x**2 / eps)) / \n",
    "                (erf(1 / k) + k / pi**0.5 * np.exp(-0.5 / eps)))    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 50,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "class Problem10(TestProblem):\n",
    "    def __init__(self, eps):\n",
    "        super(Problem10, self).__init__(-1, 1, 0, 2, eps)\n",
    "    \n",
    "    def rhs(self, x, y):\n",
    "        return -x * y[1]\n",
    "    \n",
    "    def ref(self, x):\n",
    "        k = (2 * self.eps) ** 0.5\n",
    "        return 1 + erf(x / k) / erf(1 / k)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 51,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "class Problem12(TestProblem):\n",
    "    def __init__(self, eps):\n",
    "        super(Problem12, self).__init__(-1, 1, -1, 0, eps)    \n",
    "    \n",
    "    def rhs(self, x, y):\n",
    "        eps = self.eps\n",
    "        return y[0] - (1 + eps * pi**2) * np.cos(pi * x)\n",
    "        \n",
    "    def ref(self, x):\n",
    "        k = self.eps ** -0.5\n",
    "        return (np.cos(pi * x) + \n",
    "                (np.exp(k * (1 + x)) - np.exp(-k * (1 + x))) / \n",
    "                (np.exp(2 * k) - np.exp(-2 * k)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 52,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "class Problem13(TestProblem):\n",
    "    def __init__(self, eps):\n",
    "        super(Problem13, self).__init__(-1, 1, 0, -1, eps)\n",
    "    \n",
    "    def rhs(self, x, y):\n",
    "        eps = self.eps\n",
    "        return y[0] - (1 + eps * pi**2) * np.cos(pi * x)\n",
    "        \n",
    "    def ref(self, x):\n",
    "        k = self.eps ** -0.5\n",
    "        return np.cos(pi * x) + np.exp(-(x + 1) * k)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 53,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "class Problem14(TestProblem):\n",
    "    def __init__(self, eps):\n",
    "        super(Problem14, self).__init__(-1, 1, 0, 0, eps)\n",
    "    \n",
    "    def rhs(self, x, y):\n",
    "        eps = self.eps\n",
    "        return y[0] - (1 + eps * pi**2) * np.cos(pi * x)\n",
    "        \n",
    "    def ref(self, x):\n",
    "        k = self.eps ** -0.5\n",
    "        return np.cos(pi * x) + np.exp(-(x + 1) * k) + np.exp((x - 1) * k)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 54,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "class Problem15(TestProblem):\n",
    "    def __init__(self, eps):\n",
    "        bc_a = airy(-eps**(-1/3))[0]\n",
    "        bc_b = airy( eps**(-1/3))[0]\n",
    "        super(Problem15, self).__init__(-1, 1, bc_a, bc_b, eps)\n",
    "        \n",
    "    def rhs(self, x, y):\n",
    "        return x * y[0]\n",
    "    \n",
    "    def ref(self, x):\n",
    "        return airy(x * self.eps**(-1/3))[0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 55,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "class Problem20(TestProblem):\n",
    "    def __init__(self, eps):\n",
    "        bc_a = 1 + eps * np.log(np.cosh(-0.745 / eps))\n",
    "        bc_b = 1 + eps * np.log(np.cosh(0.255 / eps))\n",
    "        super(Problem20, self).__init__(0, 1, bc_a, bc_b, eps)\n",
    "    \n",
    "    def rhs(self, x, y):\n",
    "        return 1 - y[1]**2\n",
    "    \n",
    "    def ref(self, x):\n",
    "        eps = self.eps\n",
    "        return 1 + eps * np.log(np.cosh((x - 0.745) / eps))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Running problems"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 61,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "problems = [\n",
    "    Problem4(1e-1), Problem4(1e-2), Problem4(1e-3), Problem4(1e-4), Problem4(1e-5),\n",
    "    Problem6(1e-1), Problem6(1e-2), Problem6(1e-3), Problem6(1e-4), Problem6(1e-5),\n",
    "    Problem7(1e-1), Problem7(1e-2), Problem7(1e-3), Problem7(1e-4), Problem7(1e-5),\n",
    "    Problem10(1e-1), Problem10(1e-2), Problem10(1e-3), Problem10(1e-4), Problem10(1e-5),\n",
    "    Problem12(1e-1), Problem12(1e-2), Problem12(1e-3), Problem12(1e-4), Problem12(1e-5),\n",
    "    Problem13(1e-1), Problem13(1e-2), Problem13(1e-3), Problem13(1e-4), Problem13(1e-5),\n",
    "    Problem14(1e-1), Problem14(1e-2), Problem14(1e-3), Problem14(1e-4), Problem14(1e-5),\n",
    "    Problem15(1e-1), Problem15(1e-2), Problem15(1e-3), Problem15(1e-4), Problem15(1e-5),\n",
    "    Problem20(1e-1), Problem20(1e-2)\n",
    "]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 62,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "     name            eps          n_nodes        success         error     \n",
      "   Problem4         1e-01           23              1          6.65e-06    \n",
      "   Problem4         1e-02           67              1          2.11e-06    \n",
      "   Problem4         1e-03           214             1          2.62e-07    \n",
      "   Problem4         1e-04           722             1          1.06e-07    \n",
      "   Problem4         1e-05           973             0          1.71e-02    \n",
      "---------------------------------------------------------------------------\n",
      "   Problem6         1e-01           22              1          5.47e-06    \n",
      "   Problem6         1e-02           51              1          2.22e-06    \n",
      "   Problem6         1e-03           108             1          3.75e-07    \n",
      "   Problem6         1e-04           276             1          1.38e-07    \n",
      "   Problem6         1e-05           830             1          5.86e-08    \n",
      "---------------------------------------------------------------------------\n",
      "   Problem7         1e-01           21              1          3.89e-06    \n",
      "   Problem7         1e-02           35              1          6.27e-07    \n",
      "   Problem7         1e-03           77              1          1.35e-07    \n",
      "   Problem7         1e-04           163             1          4.32e-09    \n",
      "   Problem7         1e-05           463             1          1.97e-10    \n",
      "---------------------------------------------------------------------------\n",
      "   Problem10        1e-01           21              1          2.79e-06    \n",
      "   Problem10        1e-02           55              1          6.18e-07    \n",
      "   Problem10        1e-03           123             1          2.79e-07    \n",
      "   Problem10        1e-04           305             1          1.25e-07    \n",
      "   Problem10        1e-05           925             1          5.26e-08    \n",
      "---------------------------------------------------------------------------\n",
      "   Problem12        1e-01           20              1          9.04e-06    \n",
      "   Problem12        1e-02           28              1          1.17e-05    \n",
      "   Problem12        1e-03           40              1          5.09e-06    \n",
      "   Problem12        1e-04           62              1          2.66e-06    \n",
      "   Problem12        1e-05           121             1          1.35e-06    \n",
      "---------------------------------------------------------------------------\n",
      "   Problem13        1e-01           20              1          3.05e-04    \n",
      "   Problem13        1e-02           28              1          1.17e-05    \n",
      "   Problem13        1e-03           40              1          5.09e-06    \n",
      "   Problem13        1e-04           62              1          2.66e-06    \n",
      "   Problem13        1e-05           121             1          1.35e-06    \n",
      "---------------------------------------------------------------------------\n",
      "   Problem14        1e-01           19              1          5.98e-04    \n",
      "   Problem14        1e-02           33              1          1.41e-05    \n",
      "   Problem14        1e-03           55              1          5.39e-06    \n",
      "   Problem14        1e-04           87              1          3.40e-06    \n",
      "   Problem14        1e-05           169             1          1.91e-06    \n",
      "---------------------------------------------------------------------------\n",
      "   Problem15        1e-01           15              1          7.78e-06    \n",
      "   Problem15        1e-02           38              1          2.04e-05    \n",
      "   Problem15        1e-03           104             1          4.27e-04    \n",
      "   Problem15        1e-04           310             1          7.38e-05    \n",
      "   Problem15        1e-05           956             1          3.32e-04    \n",
      "---------------------------------------------------------------------------\n",
      "   Problem20        1e-01           25              1          5.48e-07    \n",
      "   Problem20        1e-02           325             1          1.62e+142   \n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/nmayorov/anaconda/lib/python3.4/site-packages/ipykernel/__main__.py:8: RuntimeWarning: overflow encountered in square\n",
      "/Users/nmayorov/PythonDev/scipy/scipy/integrate/_bvp.py:515: RuntimeWarning: invalid value encountered in true_divide\n",
      "  res_middle /= 1 + np.abs(f_middle)\n",
      "/Users/nmayorov/PythonDev/scipy/scipy/integrate/_bvp.py:516: RuntimeWarning: invalid value encountered in true_divide\n",
      "  res_1 /= 1 + np.abs(f1)\n",
      "/Users/nmayorov/PythonDev/scipy/scipy/integrate/_bvp.py:517: RuntimeWarning: invalid value encountered in true_divide\n",
      "  res_2 /= 1 + np.abs(f2)\n",
      "/Users/nmayorov/PythonDev/scipy/scipy/integrate/_bvp.py:1030: RuntimeWarning: invalid value encountered in greater\n",
      "  insert_1, = np.nonzero((rms_res > tol) &\n",
      "/Users/nmayorov/PythonDev/scipy/scipy/integrate/_bvp.py:1031: RuntimeWarning: invalid value encountered in less\n",
      "  (rms_res < 100 * tol))\n",
      "/Users/nmayorov/PythonDev/scipy/scipy/integrate/_bvp.py:1032: RuntimeWarning: invalid value encountered in greater_equal\n",
      "  insert_2, = np.nonzero(rms_res >= 100 * tol)\n"
     ]
    }
   ],
   "source": [
    "header = \"{:^15}{:^15}{:^15}{:^15}{:^15}\".format(\"name\", \"eps\", \"n_nodes\", \"success\", \"error\")\n",
    "print(header)\n",
    "name_prev = None\n",
    "for p in problems:    \n",
    "    name = p.__class__.__name__\n",
    "    if name_prev is not None and name != name_prev:\n",
    "        print(\"-\" * len(header))\n",
    "    name_prev = name\n",
    "    \n",
    "    eps = p.eps\n",
    "    success, mesh_size, error = p.run()\n",
    "    print(\"{:^15}{:^15.0e}{:^15}{:^15}{:^15.2e}\".format(name, eps, mesh_size, success, error))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.4.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}