Threshold factor calculations.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "c29b3369",
   "metadata": {},
   "source": [
    "This notebook is a demonstration of the characterization of the threshold factor (or linear SSP coefficient).  Here we test the formulas on the classical RK4 method, which has threshold factor 1.  See https://www.overleaf.com/read/vmzhftcpnmth."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "id": "3dc81301",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from nodepy import rk\n",
    "from math import comb"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "id": "a789a36d",
   "metadata": {},
   "outputs": [],
   "source": [
    "def M_mat(r,s):\n",
    "    M = np.zeros((s+1,s+1))\n",
    "    for i in range(0,s+1):\n",
    "        for j in range(i,s+1):\n",
    "            M[i,j] = comb(j,i)*r**(j-i)\n",
    "    return M"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "id": "c8dc8df5",
   "metadata": {},
   "outputs": [],
   "source": [
    "rkm = rk.loadRKM('RK44')\n",
    "s = len(rkm)\n",
    "e = np.ones(s)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "id": "dcd31436",
   "metadata": {},
   "outputs": [],
   "source": [
    "gamma = np.zeros(s+1)\n",
    "gamma[0]=1\n",
    "for i in range(1,s+1):\n",
    "    Apow = np.linalg.matrix_power(rkm.A,i-1)\n",
    "    gamma[i] = np.dot(rkm.b,np.dot(Apow,e))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "id": "13d901a5",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.0001666666666666668\n"
     ]
    }
   ],
   "source": [
    "r = 0.999\n",
    "M = M_mat(r,s)\n",
    "mu = np.linalg.solve(M,gamma)\n",
    "print(np.min(mu))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "id": "ef5e737e",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.0\n"
     ]
    }
   ],
   "source": [
    "r = 1.00\n",
    "M = M_mat(r,s)\n",
    "mu = np.linalg.solve(M,gamma)\n",
    "print(np.min(mu))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "id": "20ec8e1c",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "-0.0001666666666666483\n"
     ]
    }
   ],
   "source": [
    "r = 1.001\n",
    "M = M_mat(r,s)\n",
    "mu = np.linalg.solve(M,gamma)\n",
    "print(np.min(mu))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7acc42b1",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.10.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}