backward integration bug

SciML · Dec 23, 2023 · 0bd031e · 0bd031e
1 parent 38f1ca4
commit 0bd031e
Show file tree

Hide file tree

Showing 8 changed files with 1,064 additions and 328 deletions.
diff --git a/Examples/.ipynb_checkpoints/Simple Pendulum- Forward and Back integration-checkpoint.ipynb b/Examples/.ipynb_checkpoints/Simple Pendulum- Forward and Back integration-checkpoint.ipynb
@@ -0,0 +1,334 @@
+{
+ "cells": [
+  {
+   "cell_type": "markdown",
+   "id": "c5b2c9a2",
+   "metadata": {},
+   "source": [
+    "# Simple Pendulum example: forward and backward integrations"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "b343f609",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "using OrdinaryDiffEq \n",
+    "using IRKGaussLegendre\n",
+    "using Plots"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "fb9278c5",
+   "metadata": {},
+   "source": [
+    "## ODE definition "
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "df7272c0",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "#Constants\n",
+    "const g = 9.81\n",
+    "L = 1.0\n",
+    "\n",
+    "#Initial Conditions\n",
+    "u₀ = [0, π / 2]\n",
+    "\n",
+    "#Define the problem\n",
+    "function simplependulum(du, u, p, t)\n",
+    "    θ = u[1]\n",
+    "    dθ = u[2]\n",
+    "    du[1] = dθ\n",
+    "    du[2] = -(g / L) * sin(θ)\n",
+    "end"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "20c50ccf",
+   "metadata": {},
+   "source": [
+    "# Forward Integrations "
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "c9b43ebd",
+   "metadata": {},
+   "source": [
+    "### Case 1 (adaptive)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "028c6101",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "tspan = (0.0, 6.3)\n",
+    "\n",
+    "#Pass to solvers\n",
+    "prob = ODEProblem(simplependulum, u₀, tspan)\n",
+    "#sol = solve(prob, Tsit5())\n",
+    "sol = solve(prob, IRKGL16(), reltol=1e-14, abstol=1e-14)\n",
+    "\n",
+    "#Plot\n",
+    "plot(sol, linewidth = 2, title = \"Simple Pendulum Problem\", xaxis = \"Time\",\n",
+    "    yaxis = \"Height\", label = [\"\\\\theta\" \"d\\\\theta\"])"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "4ffb595f",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "@show (sol.retcode, sol.t[1], sol.t[end]);\n",
+    "sol.t"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "bed51f6d",
+   "metadata": {},
+   "source": [
+    "### Case 2 (adaptive)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "3c6f59bb",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "tspan = (-1., 6.3)\n",
+    "\n",
+    "#Pass to solvers\n",
+    "prob = ODEProblem(simplependulum, u₀, tspan)\n",
+    "#sol = solve(prob, Tsit5())\n",
+    "sol = solve(prob, IRKGL16(), reltol=1e-14, abstol=1e-14)\n",
+    "\n",
+    "#Plot\n",
+    "plot(sol, linewidth = 2, title = \"Simple Pendulum Problem\", xaxis = \"Time\",\n",
+    "    yaxis = \"Height\", label = [\"\\\\theta\" \"d\\\\theta\"])"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "507ae9bf",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "@show (sol.retcode, sol.t[1], sol.t[end]);\n",
+    "sol.t"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "33e6acfa",
+   "metadata": {},
+   "source": [
+    "### Case 3 (constant step size)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "596eb99e",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "tspan = (-1., 6.3)\n",
+    "\n",
+    "dt0=0.1\n",
+    "\n",
+    "#Pass to solvers\n",
+    "prob = ODEProblem(simplependulum, u₀, tspan)\n",
+    "#sol = solve(prob, Tsit5())\n",
+    "sol = solve(prob, IRKGL16(), dt=dt0, adaptive=false)\n",
+    "\n",
+    "#Plot\n",
+    "plot(sol, linewidth = 2, title = \"Simple Pendulum Problem\", xaxis = \"Time\",\n",
+    "    yaxis = \"Height\", label = [\"\\\\theta\" \"d\\\\theta\"])"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "688c60a3",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "@show (sol.retcode, sol.t[1], sol.t[end]);\n",
+    "sol.t"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "cf87f1ca",
+   "metadata": {},
+   "source": [
+    "# Backward Integrations"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "fb08bb81",
+   "metadata": {},
+   "source": [
+    "### Case 1"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "bc2a51d9",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "tspan = (6.3, 0.)\n",
+    "\n",
+    "#Pass to solvers\n",
+    "prob = ODEProblem(simplependulum, u₀, tspan)\n",
+    "sol = solve(prob, IRKGL16(), reltol=1e-14, abstol=1e-14)\n",
+    "\n",
+    "#Plot\n",
+    "plot(sol, linewidth = 2, title = \"Simple Pendulum Problem\", xaxis = \"Time\",\n",
+    "    yaxis = \"Height\", label = [\"\\\\theta\" \"d\\\\theta\"])"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "05c59aae",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "@show (sol.retcode, sol.t[1], sol.t[end]);\n",
+    "sol.t"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "c5a72290",
+   "metadata": {},
+   "source": [
+    "**Note**\n",
+    "\n",
+    "-When tf=0, the integrator not return exactly zero\n",
+    "\n",
+    "-If tf=1., the integrator return exactly one"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "a17c971e",
+   "metadata": {},
+   "source": [
+    "### Case 2"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "e050fa3d",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "tspan = (0.0, -6.3)\n",
+    "\n",
+    "#Pass to solvers\n",
+    "prob = ODEProblem(simplependulum, u₀, tspan)\n",
+    "sol = solve(prob, IRKGL16(), reltol=1e-14, abstol=1e-14)\n",
+    "\n",
+    "#Plot\n",
+    "plot(sol, linewidth = 2, title = \"Simple Pendulum Problem\", xaxis = \"Time\",\n",
+    "    yaxis = \"Height\", label = [\"\\\\theta\" \"d\\\\theta\"])"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "cf5ae474",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "@show (sol.retcode, sol.t[1], sol.t[end]);\n",
+    "sol.t"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "bf0b8d02",
+   "metadata": {},
+   "source": [
+    "### Case 3"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "840c2638",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "tspan = (6.3, 0.)\n",
+    "\n",
+    "dt0=0.1\n",
+    "#Pass to solvers\n",
+    "prob = ODEProblem(simplependulum, u₀, tspan)\n",
+    "sol = solve(prob, IRKGL16(), dt=dt0, adaptive=false)\n",
+    "\n",
+    "#Plot\n",
+    "plot(sol, linewidth = 2, title = \"Simple Pendulum Problem\", xaxis = \"Time\",\n",
+    "    yaxis = \"Height\", label = [\"\\\\theta\" \"d\\\\theta\"])"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "bc9fefd2",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "@show (sol.retcode, sol.t[1], sol.t[end]);\n",
+    "sol.t"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "197b6d79",
+   "metadata": {},
+   "outputs": [],
+   "source": []
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "Julia 1.9.1",
+   "language": "julia",
+   "name": "julia-1.9"
+  },
+  "language_info": {
+   "file_extension": ".jl",
+   "mimetype": "application/julia",
+   "name": "julia",
+   "version": "1.9.1"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 5
+}