{ "cells": [ { "cell_type": "markdown", "id": "1fe536ed", "metadata": {}, "source": [ "# Solving an ODE with a forcing term" ] }, { "cell_type": "markdown", "id": "598ab169-05d8-4733-a6cc-9fa91aa92198", "metadata": {}, "source": [ "This example demonstrates how to incorporate an external forcing term into the solve. This is really simple: just evaluate it as part of the vector field like anything else.\n", "\n", "This example is available as a Jupyter notebook [here](https://github.com/patrick-kidger/diffrax/blob/main/examples/forcing.ipynb)." ] }, { "cell_type": "code", "execution_count": 1, "id": "6d6bdf63", "metadata": { "tags": [] }, "outputs": [], "source": [ "import jax\n", "import jax.numpy as jnp\n", "import matplotlib.pyplot as plt\n", "from diffrax import diffeqsolve, ODETerm, SaveAt, Tsit5\n", "\n", "\n", "def force(t, args):\n", " m, c = args\n", " return m * t + c\n", "\n", "\n", "def vector_field(t, y, args):\n", " return -y + force(t, args)\n", "\n", "\n", "@jax.jit\n", "def solve(y0, args):\n", " term = ODETerm(vector_field)\n", " solver = Tsit5()\n", " t0 = 0\n", " t1 = 10\n", " dt0 = 0.1\n", " saveat = SaveAt(ts=jnp.linspace(t0, t1, 1000))\n", " sol = diffeqsolve(term, solver, t0, t1, dt0, y0, args=args, saveat=saveat)\n", " return sol\n", "\n", "\n", "y0 = 1.0\n", "args = (0.1, 0.02)\n", "sol = solve(y0, args)" ] }, { "cell_type": "code", "execution_count": 2, "id": "9654fd84-19b9-4a0b-bff6-d20f36c4f333", "metadata": { "tags": [] }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjcAAAGwCAYAAABVdURTAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAABJhUlEQVR4nO3de1yT5/3/8XcSIOEYlEM4COKxSlVAEES7HlY61zrXdlur1qrFiq2zWy2/Hepa9Tu7ytZ9a7u1tirV1trWQw/reXaOb0+2CIpitZ4FBYUEEEkgQIDk/v0RiFJRUUmu5M77+XjkD2MiH3iIvLzu675vhSRJEoiIiIhkQil6ACIiIqK+xLghIiIiWWHcEBERkawwboiIiEhWGDdEREQkK4wbIiIikhXGDREREcmKj+gBXM1ms6GqqgrBwcFQKBSixyEiIqJekCQJjY2NiImJgVJ56bUZr4ubqqoqxMXFiR6DiIiIrkJlZSUGDBhwydd4XdwEBwcDsH9xQkJCBE9DREREvWEymRAXF+f4OX4pXhc3XYeiQkJCGDdEREQepjdbSrihmIiIiGSFcUNERESywrghIiIiWWHcEBERkawwboiIiEhWGDdEREQkK4wbIiIikhXGDREREckK44aIiIhkhXFDREREsiI0br766itMmTIFMTExUCgUeP/99y/7ni+++AJjx46FWq3G0KFD8dprrzl9TiIiIvIcQuPGbDYjKSkJK1eu7NXry8vLMXnyZNxyyy0oLS3FwoULMXfuXHz22WdOnpSIiIg8hdAbZ95+++24/fbbe/36VatWYdCgQXj22WcBACNHjsT27dvx3HPPYdKkSc4as1dsNgl1ZguaLVYkhAcKnYWIiMibedSem8LCQmRlZXV7btKkSSgsLLzoeywWC0wmU7eHM3xzvA7pTxdg3oZdTvnziYiIqHc8Km70ej10Ol2353Q6HUwmE1paWnp8T15eHrRareMRFxfnlNmiQjT2GY2tTvnziYiIqHc8Km6uxqJFi2A0Gh2PyspKp3wcndYeN6bWDrS0WZ3yMYiIiOjyhO65uVJRUVEwGAzdnjMYDAgJCYG/v3+P71Gr1VCr1U6fLVjtgwA/FZrbrNCbWjGI+26IiIiE8KiVm8zMTBQUFHR7btu2bcjMzBQ00TkKhYKHpoiIiNyA0LhpampCaWkpSktLAdhP9S4tLUVFRQUA+yGlWbNmOV7/8MMPo6ysDH/4wx9w6NAhvPTSS9iyZQsee+wxEeNfQNcZNwYT44aIiEgUoXGza9cupKSkICUlBQCQm5uLlJQULFmyBABQXV3tCB0AGDRoED755BNs27YNSUlJePbZZ/HKK68IPw28S1Tnvhs944aIiEgYoXtubr75ZkiSdNHf7+nqwzfffDP27NnjxKmuno6HpYiIiITzqD037i4qxL5xmYeliIiIxGHc9CEeliIiIhKPcdOHHBuKeViKiIhIGMZNH+paualptMBmu/heIiIiInIexk0fighSQ6kAOjpvoklERESux7jpQz4qJcKDOjcVGxk3REREIjBu+hg3FRMREYnFuOljjmvdMG6IiIiEYNz0sSieMUVERCQU46aP8bAUERGRWIybPsabZxIREYnFuOljUby/FBERkVCMmz4WpbWfCs7DUkRERGIwbvpY12GpxtYONLd1CJ6GiIjI+zBu+liwxheBfioAPDRFREQkAuPGCXQ8Y4qIiEgYxo0TcFMxERGROIwbJ4jiVYqJiIiEYdw4QdeF/HiVYiIiItdj3DgBr1JMREQkDuPGCc7dPNMieBIiIiLvw7hxAt48k4iISBzGjRN0HZaqbbLAapMET0NERORdGDdOEB6kho9SAatNQm0jD00RERG5EuPGCVRKhWPfzemGFsHTEBEReRfGjZNEdx6aqjYyboiIiFyJceMkMaH+AIDqBm4qJiIiciXGjZNEh9pXbqq4ckNERF7k1NlmSJLYk2kYN04So+XKDREReZf9p42Y/M/tWPrh90IDh3HjJF17brhyQ0RE3uD7KiPuX1sEY0s79p02orXdJmwWxo2TdO25qeLKDRERydz3VUbMeKUIDc3tSIkPxetz0uHvpxI2D+PGSbripq7JAkuHVfA0REREznF+2CTHhWL9nHQEa3yFzsS4cZJ+Ab5Q+9i/vAYjL+RHRETyc6DK1C1sXn8wHSGCwwZg3DiNQqE4d2iK+26IiEhm7GGzAw3N7Uhyo7ABGDdO5dhUzKsUExGRjBystofN2eZ2JA3Q4vU57hM2gBvEzcqVK5GQkACNRoOMjAwUFxdf9LXt7e1YtmwZhgwZAo1Gg6SkJGzdutWF014Zx4X8eHdwIiKSiUN6+6EoR9g8mAGtv/uEDSA4bjZv3ozc3FwsXboUu3fvRlJSEiZNmoSampoeX//kk09i9erVeOGFF3DgwAE8/PDDuPvuu7Fnzx4XT947MVy5ISIiGTmkN+G+/CLUm9swxk3DBhAcNytWrEBOTg6ys7ORmJiIVatWISAgAOvWrevx9Rs2bMCf/vQn3HHHHRg8eDDmz5+PO+64A88+++xFP4bFYoHJZOr2cJVortwQEZFMHNY3dgubDW4aNoDAuGlra0NJSQmysrLODaNUIisrC4WFhT2+x2KxQKPRdHvO398f27dvv+jHycvLg1ardTzi4uL65hPoBe65ISIiObCHzQ7Um9swOlaLDXPcN2wAgXFTV1cHq9UKnU7X7XmdTge9Xt/jeyZNmoQVK1bg6NGjsNls2LZtG9577z1UV1df9OMsWrQIRqPR8aisrOzTz+NSzl3Ij3FDRESe6YjBHjZnOsPmjQczoA1w37AB3GBD8ZX4xz/+gWHDhmHEiBHw8/PDI488guzsbCiVF/801Go1QkJCuj1cpWvlxtTaAbOlw2Ufl4iIqC8cMTRi+hp72IyKDfGIsAEExk14eDhUKhUMBkO35w0GA6Kionp8T0REBN5//32YzWacPHkShw4dQlBQEAYPHuyKka9YsMYXwRofAEA1r3VDREQe5Oh5KzbXx3hO2AAC48bPzw+pqakoKChwPGez2VBQUIDMzMxLvlej0SA2NhYdHR149913ceeddzp73KvWdXdw3mOKiIg8xVFDI6bn70Bdkz1s3pybgdAAP9Fj9ZrQw1K5ubnIz8/H+vXrcfDgQcyfPx9msxnZ2dkAgFmzZmHRokWO1xcVFeG9995DWVkZvv76a/z0pz+FzWbDH/7wB1GfwmVFh9oPTXHlhoiIPMGxmkZMzy9CXVMbEqM9L2wAwEfkB586dSpqa2uxZMkS6PV6JCcnY+vWrY5NxhUVFd3207S2tuLJJ59EWVkZgoKCcMcdd2DDhg0IDQ0V9BlcXtem4tNcuSEiIjd3rKYR09YUoa7J4rFhAwAKSZIk0UO4kslkglarhdFodMnm4hf/7yj+9z9HcE/qAPz9niSnfzwiIqKrcaymCdPW7OgWNv0C3SdsruTnt0edLeWJorW8kB8REbm3YzVNnXtsLBjphmFzpRg3Tta154Z3BiciInd0vNYeNrWNFoyICvb4sAEYN07XdbZUdUMrvOwIIBERubnjtfZDUV1h81bOePT38LABGDdOF9V5Ib+WdisamtsFT0NERGR3vLYJ02UYNgDjxuk0viqEB9n/svDQFBERuYOyzrCpOe9QlFzCBmDcuERs1+ngZxk3REQkVnmdGdPz7WFznc4eNmFBatFj9SnGjQsM6BcAADjFuCEiIoHK68yYtqYQBpMFw3VBeDNHfmEDMG5cYkA/+8oN44aIiEQ5UWfG9DU7HGHzVs54hMswbADGjUvE9uu6SnGz4EmIiMgbnagzY9qaHdCbWmUfNgDjxiW4ckNERKKc6Nxjoze1Ylik/MMGYNy4RGwo99wQEZHrnTxjD5tqo/eEDcC4cYmuw1LGlnY0tvJaN0RE5Hwnz9gPRVUbWzG0M2wiguUfNgDjxiWC1D4IDfAFAJxu4OoNERE5V8WZZkzvDJshEYF4KyfDa8IGYNy4TNe+G17rhoiInKniTDOmrSlEVWfYbJw3HpHBGtFjuRTjxkUGcN8NERE5WWV9M6bn70CVsRWDIwKxMcf7wgZg3LhMrOOMKZ4OTkREfa+yvhnT1uzA6YYWDI4IxKac8YgM8b6wARg3LuM4LMU9N0RE1Me6hU24d4cNwLhxma77S/GwFBER9aULwmaed4cNwLhxGd5fioiI+hrDpmeMGxfp2nNTb25Dc1uH4GmIiMjT/TBsNjJsHBg3LqL190WwxgcATwcnIqJr01PY6Bg2DowbF3IcmuKmYiIiukrnh80ghk2PGDcuxE3FRER0LbquY9MVNpsYNj1i3LjQAF7rhoiIrtKps/awOXW2c8Umh2FzMYwbF+ItGIiI6GqcOms/FHV+2ERpGTYXw7hxoXMrN4wbIiLqnfPDJiEsgGHTC4wbF+raUMyrFBMRUW+cfygqISwAm+ZlMmx6gXHjQl0bimsbLWhttwqehoiI3NnphhZMz9+ByvrOFZt5XLHpLcaNC4UG+CLQTwWAqzdERHRxpxtaMG1NISrrWzCwM2yitf6ix/IYjBsXUigUvA0DERFd0g/DZhPD5ooxblyMp4MTEdHFnG5owfQ1Oxg214hx42Jx/e0rNxX1jBsiIjqnqjNsKuqb7Yeichg2V4tx42LxXXFzhnFDRER2VQ0tmNYZNvH97WETE8qwuVqMGxcbGMaVGyIiOueHYbNpHsPmWgmPm5UrVyIhIQEajQYZGRkoLi6+5Ouff/55XHfddfD390dcXBwee+wxtLa2umjaa3f+yo0kSYKnISIikRg2ziE0bjZv3ozc3FwsXboUu3fvRlJSEiZNmoSampoeX//WW2/h8ccfx9KlS3Hw4EGsXbsWmzdvxp/+9CcXT371us6WarR0wNjSLngaIiISparzOjYMm74nNG5WrFiBnJwcZGdnIzExEatWrUJAQADWrVvX4+u//fZbTJw4Effddx8SEhLwk5/8BNOnT7/kao/FYoHJZOr2EMnfT4XIYDUAHpoiIvJW1UZ72Jw807nHhmHTp4TFTVtbG0pKSpCVlXVuGKUSWVlZKCws7PE9EyZMQElJiSNmysrK8Omnn+KOO+646MfJy8uDVqt1POLi4vr2E7kKXYemTnJTMRGR16k22g9FnTzTjLj+/tg4b7zjCvbUN4TFTV1dHaxWK3Q6XbfndTod9Hp9j++57777sGzZMtxwww3w9fXFkCFDcPPNN1/ysNSiRYtgNBodj8rKyj79PK5GPDcVExF5pR+GzaZ5mQwbJxC+ofhKfPHFF1i+fDleeukl7N69G++99x4++eQTPPXUUxd9j1qtRkhISLeHaF0rN5WMGyIir1FttF/HhmHjfD6iPnB4eDhUKhUMBkO35w0GA6Kionp8z+LFizFz5kzMnTsXADB69GiYzWbMmzcPTzzxBJRKz2i1eF7Ij4jIq+iNrZi+ZgdOdB2KyuGhKGcSVgN+fn5ITU1FQUGB4zmbzYaCggJkZmb2+J7m5uYLAkalst+I0pNOq2bcEBF5D72xFdPWFOLEmWYM6GcPm64zZ8k5hK3cAEBubi5mz56NtLQ0pKen4/nnn4fZbEZ2djYAYNasWYiNjUVeXh4AYMqUKVixYgVSUlKQkZGBY8eOYfHixZgyZYojcjxBV9xUNbSg3WqDr8ozVpyIiOjK/DBsNs1j2LiC0LiZOnUqamtrsWTJEuj1eiQnJ2Pr1q2OTcYVFRXdVmqefPJJKBQKPPnkkzh9+jQiIiIwZcoUPP3006I+hasSEayGxleJ1nYbTp9tQUJ4oOiRiIioj+mNrZiev4NhI4BC8qTjOX3AZDJBq9XCaDQK3Vz8k+e+xBFDE16fk44bh0cIm4OIiPpeV9iU15kdh6K6bpxMV+dKfn7zeIgg3HdDRCRP54dNbCjDRgTGjSBxPB2ciEh2DKbuYbNpHsNGBMaNIFy5ISKSF4OpFdPWMGzcAeNGkIFhvAUDEZFcGEz269gwbNwD40aQ869S7GV7uomIZKUrbMoYNm6DcSNI1+mAjZYONDS3C56GiIiuRg3Dxi0xbgTR+KqgC1ED4L4bIiJPVNO5x4Zh434YNwJxUzERkWeqMbViWv65sOHp3u6FcSNQfH/7lYkZN0REnsMRNrXnwiY+jGHjThg3AnWdMVVeZxY8CRER9UZN53VsymrNiNFqGDZuinEjUNc9pU4wboiI3F5Noz1sjneGzaZ5mQwbN8W4EWhQWGfcnGHcEBG5s5pG+1lRDBvPwLgRKCHc/o1R19SGxlaeDk5E5I5+GDYb5/FQlLtj3AgUrPFFeJAfAOBEHTcVExG5m/PDJrozbAZ2rrqT+2LcCJbQ+U1SzkNTRERupaaxFfflFznCZhPDxmMwbgTjpmIiIvdT22jBfflFOFbTxLDxQIwbwQYxboiI3EptowXT83c4wmZjDsPG0zBuBONhKSIi93F+2ESF2MOma4WdPAfjRrCuM6a4ckNEJJb9UNS5sNk0j2HjqRg3gnWt3JxtboeRdwcnIhKiK2yOMmxkgXEjWKDaB5HB9ruD89AUEZHr1TV1D5uNDBuPx7hxA13fROV1TYInISLyLnVNFkxfYw8bXYgaG+eNd5zoQZ6LceMGum7DUM4L+RERucz5Kza6EDU2zctk2MgE48YNDIrg6eBERK7UFTZHDAwbOWLcuIEE3kCTiMhl6posmJFf5AibjTk8FCU3jBs3MMix58YMSZIET0NEJF9dYXPY0IjIYHvYDI4IEj0W9THGjRsY2Hl32cbWDtSb2wRPQ0QkT2d+EDab5jFs5Ipx4wY0virEaDUAeGiKiMgZzjTZ7xXFsPEOjBs3ce50cJ4xRUTUl34YNhsZNrLHuHETvDs4EVHfO9NkwYxXuofNEIaN7DFu3MS5a90wboiI+kJX2BzSM2y8DePGTQyJtMfN8VpepZiI6FrVm9scYRPBsPE6jBs30fVNV15nhtXG08GJiK5WvbkN9+XvcITNJoaN12HcuIkB/QLgp1LC0mHD6bMtoschIvJIPwybjTkMG2/kFnGzcuVKJCQkQKPRICMjA8XFxRd97c033wyFQnHBY/LkyS6cuO+plArHxfx4aIqI6Mr1FDZDIxk23kh43GzevBm5ublYunQpdu/ejaSkJEyaNAk1NTU9vv69995DdXW147F//36oVCrcc889Lp6873HfDRHR1Tl/j014EMPG2wmPmxUrViAnJwfZ2dlITEzEqlWrEBAQgHXr1vX4+v79+yMqKsrx2LZtGwICAmQRN0M7l04ZN0REvdcVNgerTQgPsu+xYdh4N6Fx09bWhpKSEmRlZTmeUyqVyMrKQmFhYa/+jLVr12LatGkIDOz5pmcWiwUmk6nbw10N6fxmPF7D08GJiHrj7AVhk8GwIbFxU1dXB6vVCp1O1+15nU4HvV5/2fcXFxdj//79mDt37kVfk5eXB61W63jExcVd89zOMoQrN0REvXbW3Ib7LgibYNFjkRsQfljqWqxduxajR49Genr6RV+zaNEiGI1Gx6OystKFE16Zrg3FZ8xtOMsbaBIRXdQPV2w25jBs6ByhcRMeHg6VSgWDwdDteYPBgKioqEu+12w2Y9OmTXjwwQcv+Tq1Wo2QkJBuD3cVqPZx3ECzrI6rN0REPanvXLE5UG1CeJAfNuZkYJiOYUPnCI0bPz8/pKamoqCgwPGczWZDQUEBMjMzL/net99+GxaLBffff7+zx3Sprn03x2oYN0REP9R1uvf5m4cZNvRDwg9L5ebmIj8/H+vXr8fBgwcxf/58mM1mZGdnAwBmzZqFRYsWXfC+tWvX4q677kJYWJirR3aqc/tuuKmYiOh851/Hhnts6FJ8RA8wdepU1NbWYsmSJdDr9UhOTsbWrVsdm4wrKiqgVHZvsMOHD2P79u34z3/+I2Jkpzp3xhRXboiIupx/E0ye7k2Xo5AkyatuZGQymaDVamE0Gt1y/823x+twX34REsIC8MXvbxE9DhGRcOeHDa887L2u5Oe38MNS1F3Xhfwq6pth6bAKnoaISKwzTRbcl8+woSvDuHEzEcFqBKt9YJOAk2eaRY9DRCRMXWfYHDY0IjKYh6Ko9xg3bkahUHDfDRF5PXvY7HCEzcZ5vLs39R7jxg11fQPzdHAi8kZdYXPE0ARdiH3FhmFDV0L42VJ0Id4dnIi8VW2jPWyO1nSFTabj6u1EvcWVGzfkWLlh3BCRFzk/bKJCNAwbumqMGzc07LyrFNtsXnWmPhF5qZrGVkw/L2w2zhvPsKGrxrhxQwPDAuHno0Rruw2VZ3nGFBHJW01jK6av2YFjjhUbhg1dG8aNG1IpFY7r3RzWNwqehojIebrC5nitGdFae9gkMGzoGjFu3NR1Ufb7pRzlGVNEJFM1JoYNOQfjxk0N77zLLVduiEiOakytmJZvD5uYzrAZGMawob7BU8Hd1HCd/bDUEQPjhojkpStsyhxhk4n4sADRY5GMcOXGTXWt3JTVmtFhtQmehoiobxhMrZi2xh42saH+DBtyCsaNm4oN9UeAnwptVhtO8B5TRCQDhs49NmV1XWEznmFDTsG4cVNKpQLDOldveGiKiDyd3ti5YnNe2MT1Z9iQczBu3NjwSO67ISLPpzfaL9BXzrAhF+GGYjfWdTo444aIPJV9xaYQJ840M2zIZbhy48aG8XRwIvJg1cYWR9gM6MewIddh3Lix6zrj5sSZZlg6rIKnISLqPXvY7GDYkBCMGzemC1EjROMDq01CWa1Z9DhERL1S1WAPm5NnmhHX3x42A/oxbMh1GDduTKFQOK53w303ROQJLgybTIYNuRzjxs0N56ZiIvIQpzvDpqK+GfH9A7BpXiZiQ/1Fj0VeiGdLublzp4PzBppE5L7sYVOIyvoWxPcPwMZ54xk2JMwVr9zMnj0bX331lTNmoR5w5YaI3N2ps83dwmYTw4YEu+K4MRqNyMrKwrBhw7B8+XKcPn3aGXNRp64zpirqm2G2dAiehoiou1NnmzE9fwcq61swMMweNjEMGxLsiuPm/fffx+nTpzF//nxs3rwZCQkJuP322/HOO++gvb3dGTN6tbAgNSKD1ZAk4DBXb4jIjdhXbBg25H6uakNxREQEcnNzsXfvXhQVFWHo0KGYOXMmYmJi8Nhjj+Ho0aN9PadXGxkdAgA4UGUSPAkRkV1lvT1sTp1tQUJn2ERrGTbkHq7pbKnq6mps27YN27Ztg0qlwh133IF9+/YhMTERzz33XF/N6PW64uZgNeOGiMS7MGwyGTbkVq44btrb2/Huu+/iZz/7GQYOHIi3334bCxcuRFVVFdavX4///ve/2LJlC5YtW+aMeb1SYgzjhojcQ1fYnG5owaDwQGyal4korUb0WETdXPGp4NHR0bDZbJg+fTqKi4uRnJx8wWtuueUWhIaG9sF4BACJ0fZNxYf0jbDZJCiVCsETEZE3+mHYbMwZz7Aht3TFcfPcc8/hnnvugUZz8b/QoaGhKC8vv6bB6JyEsECofZRobrPiZH0zBoUHih6JiLxMxRn7WVHnVmzGQxfCsCH3dMWHpWbOnHnJsKG+56NSYkTn9W54aIqIXK3ijP06NqcbWjCYYUMegLdf8BA8Y4qIROgKmypjKwaHB2Ijw4Y8AG+/4CF4xhQRudrJM2ZMX7PDHjYRgdiUMx6RDBvyAMJXblauXImEhARoNBpkZGSguLj4kq9vaGjAggULEB0dDbVajeHDh+PTTz910bTi8IwpInKlE3VmTOsMmyEMG/IwQlduNm/ejNzcXKxatQoZGRl4/vnnMWnSJBw+fBiRkZEXvL6trQ233XYbIiMj8c477yA2NhYnT570ijOzuvbcVBlb0dDchtAAP8ETEZFcldfZV2z0JnvYbGTYkIcRunKzYsUK5OTkIDs7G4mJiVi1ahUCAgKwbt26Hl+/bt061NfX4/3338fEiRORkJCAm266CUlJSS6e3PWCNb6I62+/SNYBrt4QkZMcr23C1NWF0JtaMSwyCBvnMWzI8wiLm7a2NpSUlCArK+vcMEolsrKyUFhY2ON7PvzwQ2RmZmLBggXQ6XQYNWoUli9fDqvVetGPY7FYYDKZuj08VaJj3w3vMUVEfe9YTSOmrdmBmkYLhus6wyaYYUOeR1jc1NXVwWq1QqfTdXtep9NBr9f3+J6ysjK88847sFqt+PTTT7F48WI8++yz+Mtf/nLRj5OXlwetVut4xMXF9enn4Uo8Y4qInOWooRHT1hShttGCEVHB2JgzHuFBatFjEV0V4RuKr4TNZkNkZCTWrFmD1NRUTJ06FU888QRWrVp10fcsWrQIRqPR8aisrHThxH2LZ0wRkTMc1ttXbOqaLBgZHYK3csYjjGFDHkzYhuLw8HCoVCoYDIZuzxsMBkRFRfX4nujoaPj6+kKlUjmeGzlyJPR6Pdra2uDnd+EmW7VaDbVaHt+kXYeljtU0oa3DBj8fj2pTInJDB6tNmPFKEerNbbg+JgRvzs3gCQvk8YT9dPTz80NqaioKCgocz9lsNhQUFCAzM7PH90ycOBHHjh2DzWZzPHfkyBFER0f3GDZyM6CfP4I1Pmiz2nCspkn0OETk4b6vMuK+/B2oN7dhdKwWb80dz7AhWRD6X//c3Fzk5+dj/fr1OHjwIObPnw+z2Yzs7GwAwKxZs7Bo0SLH6+fPn4/6+no8+uijOHLkCD755BMsX74cCxYsEPUpuJRCocCoGC0AYP9po+BpiMiT7T9txIxXinC2uR1JA7R4Y24GtAG+osci6hNCr3MzdepU1NbWYsmSJdDr9UhOTsbWrVsdm4wrKiqgVJ7rr7i4OHz22Wd47LHHMGbMGMTGxuLRRx/FH//4R1GfgsuNGaBFYdkZfHe6AfeO89zN0UQkzr5TRsx4ZQdMrR1IjgvF6w+mI0TDsCH5UEiSJIkewpVMJhO0Wi2MRiNCQkJEj3PFPtpbhd9s3IOkuFB8sGCi6HGIyMPsrWzA/WuL0NjagbHxoVg/Jx3BDBvyAFfy85s7Uj3M6Fj7YamD1Sa0W22XeTUR0Tl7Ks7i/lfsYZM2sB9efzCDYUOyxLjxMAPDAuybijtsOGLgxfyIqHdKTp7FzLXFaLR0ID2hP16bk44gNe+dTPLEuPEwCoXCsXqz7xQ3FRPR5e06UY9Za4vQZOlAxqD+eDV7HMOGZI1x44FGD+iMG54xRUSXUVxej1nrimFusyJzcBhezR6HQIYNyRz/hnugrpUbng5ORJeyo+wM5ry2E81tVtwwNBz5s9Lg76e6/BuJPBxXbjzQmNhQAPYbaLZ1cFMxEV3o2+N1yH7VHjY/GhaOV2YzbMh7MG48UFx/f2j9fdFm5aZiIrrQN8fqMOe1nWhpt+Km4RHIn5UGjS/DhrwH48YDKRQKjIq1n+PPfTdEdL6vjtRizms70dpuwy3XRWD1zFSGDXkdxo2HGt15aIpxQ0Rdvjhcg7mv74Klw4askZFYxbAhL8UNxR6Kp4MT0fk+P1SDhzaUoM1qw22JOqy8byz8fPj/V/JO/JvvocZ0ng5+WM9NxUTeruCgwRE2k65n2BDxb7+HGtDv3Kbiw3puKibyVv/5Xo+H37CHzR2jo/Aiw4aIceOpFAqFY/Vm76kGscMQkRBb9+vx6zd3o90qYfKYaPxjWgp8VfxnnYjfBR4sJb4fAGBPRYPYQYjI5f69rxqPvLUbHTYJP0+KwT+mJjNsiDrxO8GDpcSHAgD2VJ4VOwgRudTH31XhkY170GGTcFdyDFbcmwQfhg2RA78bPFjygFAAQFmtGQ3NbWKHISKX+KD0NH67cQ+sNgm/GBuLZ+9NZtgQ/QC/IzxYv0A/DA4PBADsqWwQOwwROd2/9pzCY5tLYZOAX6UOwN9/lQSVUiF6LCK3w7jxcMldh6a474ZI1t4tOYXcLXthk4CpaXF45pdjGDZEF8G48XDnNhVz3w2RXG0qrsDv3tkLSQKmp8cj7xejoWTYEF0U48bDpcSFAgBKKxtgs0lihyGiPrdhx0k8/t4+SBIwK3Mgnr5rFMOG6DIYNx5uRFQwNL5KNLZ2oKyuSfQ4RNSH1m0vx+L39wMA5kwchD///HqGDVEvMG48nI9KiTGdZ03t5r4bItlY89VxLPv4AADgoZsGY/HPRkKhYNgQ9QbjRgZSuKmYSFZWfn4Myz89BAD4zY+H4vGfjmDYEF0B3hVcBlLiuKmYSA4kScLz/z2KfxQcBQDk3jYcv711mOCpiDwPV25koGvl5oihEU2WDrHDENFVkSQJ//ufw46w+eNPRzBsiK4S40YGdCEaxIb6wyYB3/EmmkQeR5Ik5P37EFZ+fhwA8OTkkZh/8xDBUxF5LsaNTHRdzG/3SR6aIvIkkiThzx8dwJqvygAAy+68HnN/NFjwVESejXEjE2kD7ftudp5g3BB5CptNwpPv78dr354AACy/ezRmZSYInYlIDrihWCbGJfQHYF+5sdokXpadyM1ZbRIWvfcdtuw6BYUC+Nsvx+DetDjRYxHJAlduZGJkdAiC1D5otHTgkN4kehwiugSrTcLv396LLbtOQakAVtybxLAh6kOMG5lQKRUY23Voqrxe8DREdDEdVhsWbi7Fe3tOQ6VU4B/TUnB3ygDRYxHJCuNGRtITOuOGm4qJ3FK71YbfbNyDj/ZWwVelwMr7UjAlKUb0WESywz03MpLWue9mZ3k9JEniFU2J3Iilw4pH3tqDbQcM8FMp8dKMschK1Ikei0iWuHIjI8lxofBVKVDTaEFlfYvocYioU2u7FQ9vKLGHjY8Sa2alMmyInIhxIyMaXxVGx2oBAMUnuO+GyB20tFmR8/oufH64FhpfJdbNHoebr4sUPRaRrLlF3KxcuRIJCQnQaDTIyMhAcXHxRV/72muvQaFQdHtoNBoXTuvexg06d2iKiMRqbuvAnNd24uujdQjwU+G17HTcMCxc9FhEsic8bjZv3ozc3FwsXboUu3fvRlJSEiZNmoSampqLvickJATV1dWOx8mTJ104sXtL79p3c5JxQyRSk6UDD6zbicKyMwhS++D1OekYPzhM9FhEXkF43KxYsQI5OTnIzs5GYmIiVq1ahYCAAKxbt+6i71EoFIiKinI8dLqLH7u2WCwwmUzdHnKW2nk6eFmtGXVNFsHTEHknU2s7Zq4tQvGJegRrfPD6g+mODf9E5HxC46atrQ0lJSXIyspyPKdUKpGVlYXCwsKLvq+pqQkDBw5EXFwc7rzzTnz//fcXfW1eXh60Wq3jERcn7wtlhQb4YURUMACgqIyrN0Su1tDchpmvFGFPRQO0/r54a+54jI3vJ3osIq8iNG7q6upgtVovWHnR6XTQ6/U9vue6667DunXr8MEHH+CNN96AzWbDhAkTcOrUqR5fv2jRIhiNRsejsrKyzz8Pd5M5xL70/e3xOsGTEHmXuiYLpq3Zgb2njOgX4IuNOeMxeoBW9FhEXsfjrnOTmZmJzMxMx68nTJiAkSNHYvXq1XjqqacueL1arYZarXbliMJlDg7Dq9+cQGHZGdGjEHkNg6kVM14pwrGaJoQHqfFWTgaG64JFj0XklYSu3ISHh0OlUsFgMHR73mAwICoqqld/hq+vL1JSUnDs2DFnjOiRMgaHQamw77sxmFpFj0Mke6fONuPe1YU4VtOEaK0GWx4az7AhEkho3Pj5+SE1NRUFBQWO52w2GwoKCrqtzlyK1WrFvn37EB0d7awxPY7W3xfXx9iXwguPc/WGyJlO1JkxdfUOnDzTjLj+/tjyUCYGRwSJHovIqwk/Wyo3Nxf5+flYv349Dh48iPnz58NsNiM7OxsAMGvWLCxatMjx+mXLluE///kPysrKsHv3btx///04efIk5s6dK+pTcEsTOvfdMG6InOdYTSPuXV2I0w0tGBweiC0PZSKuf4DosYi8nvA9N1OnTkVtbS2WLFkCvV6P5ORkbN261bHJuKKiAkrluQY7e/YscnJyoNfr0a9fP6SmpuLbb79FYmKiqE/BLY0fEobVX5Xh2zJuKiZyhgNVJsxcW4Qz5jZcpwvGG3MzEBHsXfv7iNyVQpIkSfQQrmQymaDVamE0GhESEiJ6HKdpsnQg6c//gdUm4es/3ML/TRL1ob2VDZi1rhjGlnaMig3BhjkZ6BfoJ3osIlm7kp/fwg9LkXMEqX2Q1HkKKs+aIuo7O0/UY8YrRTC2tGNsfCjenDueYUPkZhg3MtZ1vZsd3HdD1Ce+PVaHWWuL0WTpwPjB/bHhwQxo/X1Fj0VEP8C4kbEJQ+w36PvmeB287OgjUZ/7/FANHnhtJ1rarbhxeARefSAdgWrh2xaJqAeMGxlLHdgPah8lDCYLjtY0iR6HyGNt3V+NeRt2oa3DhtsSdciflQp/P5XosYjoIhg3MqbxVSGj8y7EXx2pFTwNkWf6oPQ0Fry1B+1WCT8bE42XZoyF2odhQ+TOGDcyd+Mw+6Gpr47ylHCiK7VlZyUWbi6F1SbhV6kD8I9pKfBV8Z9NInfH71KZu3F4BACgqOwMWtutgqch8hyvF57AH979DpIE3D8+Hs/8cgxUSoXosYioFxg3MjcsMghRIRpYOmzYeaJe9DhEHmHNV8ex5IPvAQAP3jAIT905CkqGDZHHYNzInEKhwI3DOw9Ncd8N0SVJkoR//Pcoln96CADwyC1D8eTkkVAoGDZEnoRx4wV+NMx+aOpr7rshuihJkvDXrYfw3H+PAAB+95Ph+N2k6xg2RB6IceMFbhgaDoUCOKRvhMHUKnocIrdjs0l48v39WP1lGQDgyckj8ciPhwmeioiuFuPGC/QL9MOYWPutGHhoiqi7dqsNuVtK8WZRBRQKYPndozH3R4NFj0VE14Bx4yW6zpr6gnFD5NDabsWv39yN90ur4KNU4PmpybgvI170WER0jRg3XuKWEZEA7Cs37Vab4GmIxGtu68Dc9buw7YABfj5KrJ6ZijuTY0WPRUR9gHHjJZIHhCIs0A+NrR08JZy8nrGlHfe/UoTtx+oQ4KfCa9njcOtIneixiKiPMG68hFKpcKze/N/BGsHTEIlT12TBtDU7sLuiAVp/X7w5N8Nxk1kikgfGjRe5tStuDjFuyDtVNbTg3tWFOFhtQniQGpvmjUdKfD/RYxFRH2PceJEbhoXDV6VAWZ0ZZbW8Szh5lxN1ZtyzqhBltWbEhvrj7YczMTI6RPRYROQEjBsvEqzxRcYg+13CuXpD3uSwvhH3rC7E6YYWDAoPxJaHMzEoPFD0WETkJIwbL3PrSPuhqQLuuyEvUVrZgKlrClHbaMGIqGBseSgTsaH+osciIidi3HiZH3fuu9l5oh7GlnbB0xA5V+HxM5iRvwMNze1IiQ/F5nmZiAhWix6LiJyMceNlBoYFYmhkEDpsEr44zNUbkq/PD9XggVeLYW6zYsKQMLzxYAa0Ab6ixyIiF2DceKFJ19uv5/HZ93rBkxA5x0d7q5Dz+i5YOmzIGhmJdQ+MQ6DaR/RYROQijBsv9NProwEAnx+qRWu7VfA0RH1rU3EFfrtpDzpsEn6eFIOX70+FxlcleiwiciHGjRcaFRuC2FB/tLRb8SXvNUUy8vIXx/H4e/sgScD09Hg8NzUZvir+M0fkbfhd74UUCgV+OioKAPDZfh6aIs8nSRKWf3oQf9t6CAAw/+YhWH73KKiUCsGTEZEIjBsvdXtn3Gw7aEBbB2+kSZ6rw2rDH975Dmu+KgMAPHHHSPzxpyOgUDBsiLwV48ZLjY3vh4hgNRpbO1BYdkb0OERXpbXdivlv7sbbJaegVADP/GoMcm4cLHosIhKMceOllEqF46yprTw0RR6osbUdD7xajG0HDPDzUeLl+1Nxb1qc6LGIyA0wbrxY11lT//lejw4rD02R56hrsmB6/g7sKKtHkNoH67PTMen6KNFjEZGbYNx4sYzB/dE/0A9nzG08NEUe49TZZtyzqhD7T5sQFuiHTfPGI3NImOixiMiNMG68mK9K6dhY/GFpleBpiC7vqKERv3q5EOV15+7sPSpWK3osInIzjBsvd2dyLAD7vhte0I/c2Z6Ks7hndSH0plYMiwzCO/MzMTgiSPRYROSGGDdeLm1gP0RrNWi0dOCLw7ygH7mnr4/WYsYrRWhobkdyXCi2PJSJaC3v7E1EPXOLuFm5ciUSEhKg0WiQkZGB4uLiXr1v06ZNUCgUuOuuu5w7oIwplQpMSYoBYL8fD5G7+eS7asx5bSea26z40bBwvDk3A/0C/USPRURuTHjcbN68Gbm5uVi6dCl2796NpKQkTJo0CTU1l75j9YkTJ/C73/0OP/rRj1w0qXz9vDNu/nvQgMbWdsHTEJ2zofAEHtm4G+1WCZNHR+OV2Wm8ASYRXZbwuFmxYgVycnKQnZ2NxMRErFq1CgEBAVi3bt1F32O1WjFjxgz8+c9/xuDBvGDXtbo+JgSDIwJh6bBh2wGD6HGIIEkS/vezw1j8wfeQJOC+jHj8c3oK1D68ASYRXZ7QuGlra0NJSQmysrIczymVSmRlZaGwsPCi71u2bBkiIyPx4IMPXvZjWCwWmEymbg/qTqFQOFZv/rXntOBpyNu1W23447vf4cXPjwEAHssajqfv4n2iiKj3hMZNXV0drFYrdDpdt+d1Oh30+p6vmrt9+3asXbsW+fn5vfoYeXl50Gq1jkdcHK9g2pO7U+xnTW0/VodqY4vgachbNbd1YN7ru7Bll/12Cn/9xWg8mjWM94kioisi/LDUlWhsbMTMmTORn5+P8PDwXr1n0aJFMBqNjkdlZaWTp/RMA8MCkT6oPyQJeG83V2/I9erNbZieX4TPD9dC7aPE6plpmJYeL3osIvJAQnfmhYeHQ6VSwWDovs/DYDAgKurCS6kfP34cJ06cwJQpUxzP2Wz22wb4+Pjg8OHDGDJkSLf3qNVqqNVqJ0wvP/ekDkBxeT3e3lWJX988hP9bJpeprG/G7HXFKKszIzTAF2tnpyF1YH/RYxGRhxK6cuPn54fU1FQUFBQ4nrPZbCgoKEBmZuYFrx8xYgT27duH0tJSx+PnP/85brnlFpSWlvKQ0zW6Y3Q0AvxUOHGmGbtOnhU9DnmJ76uM+MXL36Ks86rD7zw8gWFDRNdE+DmVubm5mD17NtLS0pCeno7nn38eZrMZ2dnZAIBZs2YhNjYWeXl50Gg0GDVqVLf3h4aGAsAFz9OVC1T7YPLoaLxdcgpv76rEuAT+gCHn+vZYHeZtKEGTpQMjooKxfk46dCEa0WMRkYcTHjdTp05FbW0tlixZAr1ej+TkZGzdutWxybiiogJKpUdtDfJo96TF4e2SU/j4u2osnXI9rylCTvPR3irkbilFu1XC+MH9sXpmGrT+vqLHIiIZUEiSJIkewpVMJhO0Wi2MRiNCQkJEj+N2JEnCLf/7BU6cacYzvxqDe9N4qI/63rrt5Vj28QEAwOTR0VgxNYnXsCGiS7qSn99cEqFuFAoF7h1nD5o3d5wUPA3JjdUmYdlHBxxh88CEBLzAi/MRUR9j3NAF7k2Lg59Kib2njPjuVIPocUgmWtqs+PWbJVj3TTkA4PHbR2DplEQoeXE+IupjjBu6QHiQGrePtp+K/wZXb6gP1DVZMD1/Bz773gA/lRIvTE/BwzfxcgNE5ByMG+rRzPEDAQAf7q2CsZk306Srd7y2CXe/9A1KKxsQGuCLN3MyHHeiJyJyBsYN9Sh1YD+MiApGa7sN7+w+JXoc8lDF5fX4xUvforK+BfH9A/Du/Am8xAAROR3jhnqkUChwf+fqzZs7TsLLTqqjPvBB6Wnc/0oRjC3tSIkPxb9+PQFDIoJEj0VEXoBxQxd1V0osgtQ+KKsz48sjtaLHIQ8hSRJWfn4Mj24qRZvVhp9eH4WNOeMRFsTboBCRazBu6KKC1D6O69y88nW54GnIE3RYbfjTv/bh758dBgA8eMMgrJwxFhpfnupNRK7DuKFLyp6YAKUC2H6sDgeqTKLHITdmam3Hg+t3YWNxJRQK4H+mJGLxzxKh4qneRORijBu6pLj+Abh9dDQA4JXtZYKnIXdVcaYZv3zpW3x5pBYaXyVW35+KByYOEj0WEXkpxg1dVs6PBgOw3wvIYGoVPA25m50n6nHXS9/gaE0TdCFqbHkoEz+5Pkr0WETkxRg3dFnJcaEYl9AP7VYJr317QvQ45EbeLTmFGflFqDe3YVRsCD5YcAPGDAgVPRYReTnGDfXK3M7Vmzd2nISxhRf183Y2m4S/bT2E//f2XscZUVseykSUViN6NCIixg31zm0jdRiuC0JjawfWc/XGqzW3dWD+myV4+YvjAIAFtwzBSzPGIsDPR/BkRER2jBvqFaVSgUd+PAwAsHZ7OZosHYInIhGqjS341cuFjntEPTc1Cb+fNII3vyQit8K4oV6bPDoagyMCYWxpx+uFJ0SPQy5WWtmAO1/8BgeqTQgP8sPGeeNxd8oA0WMREV2AcUO9plIq8JsfDwVgv6hfcxtXb7zFOyWncO/qQtQ0WnCdLhjvL5iI1IH9RI9FRNQjxg1dkSljYjAwLAD15jZsKDwpehxysnarDf/z4ff43dt70dZhw22JOrwzPxMD+gWIHo2I6KIYN3RFfFRK/KZz781LXxznmVMydqbJgllrix2n/y/MGobV96ciWOMrdjAiostg3NAVuzslFsN1QTC2tGPVl8dFj0NO8H2VET9/8RsUlp1BoJ8Kq2emYmHWcG4cJiKPwLihK6ZSKvCHSSMAAOu2l0Nv5FWL5eTDvVX45cvf4nRDCxLCAvD+gomYxCsOE5EHYdzQVbl1ZCTGJfSDpcOGfxQcET0O9QGrTULevw/itxv3oLXdhpuGR+CDBTdgmC5Y9GhERFeEcUNXRaFQ4PHb7as3m3dW4oihUfBEdC3ONFnwwKvFWP2l/eaoD980BOseGAdtAPfXEJHnYdzQVUsd2B8/SdTBJgF//uh7SJIkeiS6CiUnz+JnL2zH10fr4O+rwgvTU/D47SOg4v4aIvJQjBu6Jk9OToSfjxLfHDuDrfv1osehKyBJEl79phxTVxei2tiKwRGBeH/BRExJihE9GhHRNWHc0DWJDwvAwzfab6r5l08OoqXNKngi6o0mSwd+s3EP/vzRAXTYJEweHY0PH7kB10Vxfw0ReT7GDV2z+TcPRWyoP043tOBlnhru9o4YGnHni9vx8XfV8FEqsORniXjxvhQEqXnjSyKSB8YNXTN/PxWenDwSALDqi+M4VtMkeCK6mA9KT+POF7/B8VozokI02PzQeMy5YRAUCu6vISL5YNxQn/jpqCjccl0E2qw2/PHd72C1cXOxOzFbOvD7t/fi0U2laGm3YuLQMHz82xuQOrC/6NGIiPoc44b6hEKhwNN3j0agnwolJ8/yruFu5PsqI6a8uB1vl5yCUgH89tZheH1OBsKD1KJHIyJyCsYN9ZmYUH88fof98NQzWw+jsr5Z8ETeretsqLtXfouyzsNQb+WMR+5tw3maNxHJGuOG+tSM9HikD+qPlnYrfvf2Xh6eEqTe3Iac13fhzx8dQJvVhqyRkfj00R9h/OAw0aMRETkd44b6lFKpwDO/HIMAPxWKyut5Y00Bvj1ehzv+8TX+e7AGfiol/mdKIvJnpaF/oJ/o0YiIXIJxQ30uITwQf/759QCAFduOYHfFWcETeYfWdiuWfXQA9+UXQW86d1G+BybybCgi8i5uETcrV65EQkICNBoNMjIyUFxcfNHXvvfee0hLS0NoaCgCAwORnJyMDRs2uHBa6o1fpQ7AlKQYWG0SHt20B6bWdtEjydp3pxow+Z9fY9035QCA6enx+Pg3NyAxJkTwZEREric8bjZv3ozc3FwsXboUu3fvRlJSEiZNmoSampoeX9+/f3888cQTKCwsxHfffYfs7GxkZ2fjs88+c/HkdCkKhQJ/uWsUYkP9UVnfgv+3ZS9s3H/T59qtNjz/3yO4+6VvcbzWjMhgNV7NHoe8X4xGgB8vykdE3kkhCb7bYUZGBsaNG4cXX3wRAGCz2RAXF4ff/OY3ePzxx3v1Z4wdOxaTJ0/GU089ddnXmkwmaLVaGI1GhITwf7XOVlrZgHtXF6Ktw4bHsobj0axhokeSjWM1jcjdshffnTICAH42JhpP3TkK/bi3hohk6Ep+fgtduWlra0NJSQmysrIczymVSmRlZaGwsPCy75ckCQUFBTh8+DBuvPHGHl9jsVhgMpm6Pch1kuNC8Ze7RgEAnvvvEfz3gEHwRJ6v3WrDys+P4Y5/bsd3p4zQ+vvin9NT8OJ9Yxk2REQQHDd1dXWwWq3Q6XTdntfpdNDrL36HaaPRiKCgIPj5+WHy5Ml44YUXcNttt/X42ry8PGi1WscjLi6uTz8Hurx70+IwK3MgAGDh5lIcqGJgXq29lQ2Y8sJ2/P2zw2jrsOGm4RH4bOGN+Dnv5E1E5CB8z83VCA4ORmlpKXbu3Imnn34aubm5+OKLL3p87aJFi2A0Gh2PyspK1w5LAIAnJydi/OD+aLJ04IFXi3HqLC/wdyWa2zrw1McHcPdL3+CQvhH9Anzx/NRkvJY9DlFajejxiIjcitAdh+Hh4VCpVDAYuh+qMBgMiIqKuuj7lEolhg4dCgBITk7GwYMHkZeXh5tvvvmC16rVaqjVvMy8aH4+SqyemYZ7VxXisKERs9cV452HJ/AwSi98eaQWT/xrH06dbQEA3JUcg8U/S0QYb59ARNQjoSs3fn5+SE1NRUFBgeM5m82GgoICZGZm9vrPsdlssFgszhiR+pDW3xevzRmHaK0Gx2vNeOC1nTC28BTxi6msb8ZDG3Zh9rpinDrbgthQf7yaPQ7PT0th2BARXYLwc0Vzc3Mxe/ZspKWlIT09Hc8//zzMZjOys7MBALNmzUJsbCzy8vIA2PfQpKWlYciQIbBYLPj000+xYcMGvPzyyyI/DeqlaK0/1s9Jx72rC7G3sgGz1hbh9QczoPX3FT2a22htt2LVl8fx8hfHYemwQaVUYHZmAnJ/MhxBauHfskREbk/4v5RTp05FbW0tlixZAr1ej+TkZGzdutWxybiiogJK5bkFJrPZjF//+tc4deoU/P39MWLECLzxxhuYOnWqqE+BrtBwXTDemjseM17Zgb2njJi5tgivz0lHaIB3H6KSJAnbDhiw7OMDjkNQmYPD8Oc7r8dwXbDg6YiIPIfw69y4Gq9z4z4OVpsw45Ui1JvbMDQyCK9lj8OAfgGixxJiT8VZ/PXfh1BUXg8AiNZq8MTkkZg8Opq3TiAiwpX9/GbckFCH9fbNxXpTKyKC1Xj1gXEYFasVPZbLlNeZ8ffPDuHTffZLH/j5KDH3hkF45MdDeYVhIqLzMG4ugXHjfqqNLch+dScO6Rvh76vCX385Gncmx4oey6n0xlas/PwYNhZXoMMmQaEAfjl2AHJvG46YUH/R4xERuR3GzSUwbtyTqbUdC97cja+P1gEAZmUOxBOTR0LtoxI8Wd863dCCVV8cx+adlWiz2gAAt1wXgT/ePgIjovj3kYjoYhg3l8C4cV9Wm4QV2w5j5efHAQCjY7X433uScF2U52+mrTjTjJe/PI53SirRbrV/y6Un9MfC24ZhwpBwwdMREbk/xs0lMG7c3/8dMuCxzXthbGmHr0qBR28dhoduGgJflWddUFuSJBSV12Pd9nJsO2hA13da5uAw/PbWYcgcEiZ2QCIiD8K4uQTGjWcwmFrxp/f2oeBQDQBgSEQgnvxZIm65LlLwZJfX0mbFJ/uq8eo35fj+vPto3Tg8Ao/cMhTpg/oLnI6IyDMxbi6BceM5JEnC+6WnseyjAzjbbL+S8Y3DI7AwaxjGxvcTPF13kiShtLIBW3adwsd7q9Bo6QAAaHyV+MXYAciekIBhvFYNEdFVY9xcAuPG8xhb2rHy82N49Ztyx36VCUPCMO/GwbhxWASUSjHXgZEkCYf0jdi6X49P9lXjWE2T4/cG9PPH9PR43Jcez/tnERH1AcbNJTBuPNfJM2as/PwY3tt9Gh02+1/b2FB//Cp1AO5KicWg8ECnz9DabsXuk2fx5dFafLZfjxNnzt3dXOOrxB2jovGrtAEYPyhMWHQREckR4+YSGDeer6qhBa98XY53Siphau1wPD80Mgi3joxE5uAwjB3YDyGaa79fVV2TBftOG7HvlBFF5Wew68RZWDpsjt/381HixmERuH1UFG67XtcnH5OIiC7EuLkExo18tLZb8dn3erxTcgqFx884VnMAQKEAhkQEYUhEIIZEBCEhLBD9Av3QL8AXWn9fx6qKJAFNlg6YWtrR0NKOqoYWnDzTjMr6ZhyvbUK1sfWCjxsZrMbEoeG4dWQkbrkuEoG8mSURkdMxbi6BcSNPxpZ2fHG4Bl8eqUXJybM4ed7homuhUACDwgMxOlaLlLhQ3DAsHEMigni/JyIiF7uSn9/8LyfJgtbfF3cmxzpu21DT2IqD1Y04XtOEsromVNa3oKG5DWeb22FsaYckSeiq+mC1D0L87Ss60VoN4vsHID4sEAPDAjAyOgRBXJkhIvIo/FebZCkyWIPIYA1uGh4hehQiInIxz7rkKxEREdFlMG6IiIhIVhg3REREJCuMGyIiIpIVxg0RERHJCuOGiIiIZIVxQ0RERLLCuCEiIiJZYdwQERGRrDBuiIiISFYYN0RERCQrjBsiIiKSFcYNERERyQrjhoiIiGTFR/QAriZJEgDAZDIJnoSIiIh6q+vndtfP8UvxurhpbGwEAMTFxQmehIiIiK5UY2MjtFrtJV+jkHqTQDJis9lQVVWF4OBgKBSKPv2zTSYT4uLiUFlZiZCQkD79s+kcfp1dg19n1+DX2XX4tXYNZ32dJUlCY2MjYmJioFReeleN163cKJVKDBgwwKkfIyQkhN84LsCvs2vw6+wa/Dq7Dr/WruGMr/PlVmy6cEMxERERyQrjhoiIiGSFcdOH1Go1li5dCrVaLXoUWePX2TX4dXYNfp1dh19r13CHr7PXbSgmIiIieePKDREREckK44aIiIhkhXFDREREssK4ISIiIllh3PSRlStXIiEhARqNBhkZGSguLhY9kuzk5eVh3LhxCA4ORmRkJO666y4cPnxY9Fiy9te//hUKhQILFy4UPYosnT59Gvfffz/CwsLg7++P0aNHY9euXaLHkhWr1YrFixdj0KBB8Pf3x5AhQ/DUU0/16v5EdHFfffUVpkyZgpiYGCgUCrz//vvdfl+SJCxZsgTR0dHw9/dHVlYWjh496rL5GDd9YPPmzcjNzcXSpUuxe/duJCUlYdKkSaipqRE9mqx8+eWXWLBgAXbs2IFt27ahvb0dP/nJT2A2m0WPJks7d+7E6tWrMWbMGNGjyNLZs2cxceJE+Pr64t///jcOHDiAZ599Fv369RM9mqz87W9/w8svv4wXX3wRBw8exN/+9jc888wzeOGFF0SP5tHMZjOSkpKwcuXKHn//mWeewT//+U+sWrUKRUVFCAwMxKRJk9Da2uqaASW6Zunp6dKCBQscv7ZarVJMTIyUl5cncCr5q6mpkQBIX375pehRZKexsVEaNmyYtG3bNummm26SHn30UdEjyc4f//hH6YYbbhA9huxNnjxZmjNnTrfnfvGLX0gzZswQNJH8AJD+9a9/OX5ts9mkqKgo6e9//7vjuYaGBkmtVksbN250yUxcublGbW1tKCkpQVZWluM5pVKJrKwsFBYWCpxM/oxGIwCgf//+gieRnwULFmDy5Mnd/l5T3/rwww+RlpaGe+65B5GRkUhJSUF+fr7osWRnwoQJKCgowJEjRwAAe/fuxfbt23H77bcLnky+ysvLodfru/37odVqkZGR4bKfi15348y+VldXB6vVCp1O1+15nU6HQ4cOCZpK/mw2GxYuXIiJEydi1KhRoseRlU2bNmH37t3YuXOn6FFkraysDC+//DJyc3Pxpz/9CTt37sRvf/tb+Pn5Yfbs2aLHk43HH38cJpMJI0aMgEqlgtVqxdNPP40ZM2aIHk229Ho9APT4c7Hr95yNcUMeacGCBdi/fz+2b98uehRZqaysxKOPPopt27ZBo9GIHkfWbDYb0tLSsHz5cgBASkoK9u/fj1WrVjFu+tCWLVvw5ptv4q233sL111+P0tJSLFy4EDExMfw6yxgPS12j8PBwqFQqGAyGbs8bDAZERUUJmkreHnnkEXz88cf4/PPPMWDAANHjyEpJSQlqamowduxY+Pj4wMfHB19++SX++c9/wsfHB1arVfSIshEdHY3ExMRuz40cORIVFRWCJpKn3//+93j88ccxbdo0jB49GjNnzsRjjz2GvLw80aPJVtfPPpE/Fxk318jPzw+pqakoKChwPGez2VBQUIDMzEyBk8mPJEl45JFH8K9//Qv/93//h0GDBokeSXZuvfVW7Nu3D6WlpY5HWloaZsyYgdLSUqhUKtEjysbEiRMvuJTBkSNHMHDgQEETyVNzczOUyu4/6lQqFWw2m6CJ5G/QoEGIiorq9nPRZDKhqKjIZT8XeViqD+Tm5mL27NlIS0tDeno6nn/+eZjNZmRnZ4seTVYWLFiAt956Cx988AGCg4Mdx261Wi38/f0FTycPwcHBF+xhCgwMRFhYGPc29bHHHnsMEyZMwPLly3HvvfeiuLgYa9aswZo1a0SPJitTpkzB008/jfj4eFx//fXYs2cPVqxYgTlz5ogezaM1NTXh2LFjjl+Xl5ejtLQU/fv3R3x8PBYuXIi//OUvGDZsGAYNGoTFixcjJiYGd911l2sGdMk5WV7ghRdekOLj4yU/Pz8pPT1d2rFjh+iRZAdAj49XX31V9GiyxlPBneejjz6SRo0aJanVamnEiBHSmjVrRI8kOyaTSXr00Uel+Ph4SaPRSIMHD5aeeOIJyWKxiB7No33++ec9/ns8e/ZsSZLsp4MvXrxY0ul0klqtlm699Vbp8OHDLptPIUm8TCMRERHJB/fcEBERkawwboiIiEhWGDdEREQkK4wbIiIikhXGDREREckK44aIiIhkhXFDREREssK4ISIiIllh3BAREZGsMG6ISDZuvvlmLFy4UPQYRCQY44aIiIhkhfeWIiJZeOCBB7B+/fpuz5WXlyMhIUHMQEQkDOOGiGTBaDTi9ttvx6hRo7Bs2TIAQEREBFQqleDJiMjVfEQPQETUF7RaLfz8/BAQEICoqCjR4xCRQNxzQ0RERLLCuCEiIiJZYdwQkWz4+fnBarWKHoOIBGPcEJFsJCQkoKioCCdOnEBdXR1sNpvokYhIAMYNEcnG7373O6hUKiQmJiIiIgIVFRWiRyIiAXgqOBEREckKV26IiIhIVhg3REREJCuMGyIiIpIVxg0RERHJCuOGiIiIZIVxQ0RERLLCuCEiIiJZYdwQERGRrDBuiIiISFYYN0RERCQrjBsiIiKSlf8Pf3L5o3sC/f0AAAAASUVORK5CYII=", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "plt.plot(sol.ts, sol.ys)\n", "plt.xlabel(\"t\")\n", "plt.ylabel(\"y\")\n", "plt.show()" ] }, { "cell_type": "markdown", "id": "2043029d-e78b-410d-9a96-5f6904f2ca05", "metadata": {}, "source": [ "Now let's consider a more complicated example: the forcing term is an interpolation, and what's more we would like to differentiate with respect to the values we are interpolating." ] }, { "cell_type": "code", "execution_count": 3, "id": "7a3f2ea2-0067-4999-a04c-5b445e7ab749", "metadata": { "tags": [] }, "outputs": [], "source": [ "from diffrax import backward_hermite_coefficients, CubicInterpolation\n", "\n", "\n", "def vector_field2(t, y, interp):\n", " return -y + interp.evaluate(t)\n", "\n", "\n", "@jax.jit\n", "@jax.grad\n", "def solve(points):\n", " t0 = 0\n", " t1 = 10\n", " ts = jnp.linspace(t0, t1, len(points))\n", " coeffs = backward_hermite_coefficients(ts, points)\n", " interp = CubicInterpolation(ts, coeffs)\n", " term = ODETerm(vector_field2)\n", " solver = Tsit5()\n", " dt0 = 0.1\n", " y0 = 1.0\n", " sol = diffeqsolve(term, solver, t0, t1, dt0, y0, args=interp)\n", " (y1,) = sol.ys\n", " return y1\n", "\n", "\n", "points = jnp.array([3.0, 0.5, -0.8, 1.8])\n", "grads = solve(points)" ] }, { "cell_type": "markdown", "id": "f34b4824-8420-4881-b5b7-78b2118de5e0", "metadata": {}, "source": [ "In this example, we computed the interpolation in advance (not repeatedly on each step!), and then just evaluated it inside the vector field." ] } ], "metadata": { "jupytext": { "formats": "ipynb,py:light" }, "kernelspec": { "display_name": "py39", "language": "python", "name": "py39" }, "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.9.16" } }, "nbformat": 4, "nbformat_minor": 5 }