Section 2.6 Starter code for Problem 11.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "ffb09b69",
   "metadata": {},
   "source": [
    "This is a Markdown field. Double-click or press return on it to edit. Press ctrl-return or shift-return to \"run\" it. Same with the code fields that follow: double-click or press return to edit, then ctrl-return or shift-return to run the code.\n",
    "\n",
    "Honors Differential Equations Fall 2023<br>\n",
    "Zill §2.6, problem #11. Obtain a numerical solution curve using Euler's and RK4.<br>\n",
    "First, use Euler's Method to approximate the Initial-Value Problem $y'=2(\\cos x)y, y(0)=1$<br>\n",
    "<br>\n",
    "Euler's Method uses the function $y_{n+1}=y_n+hf(x_n,y_n)$. Here, $f(x_n,y_n)=2(\\cos x)y$.<br>\n",
    "<br>\n",
    "Print the value of $h$ and corresponding values of $n$, $x_n$, $y_n$<br>\n",
    "Start with $h=0.25$. Then plot the solution on $0\\le x\\le 10$.\n",
    "\n",
    "To do: use RK4 and plot on the same graph."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b7cd4da3",
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import math\n",
    "from time import sleep\n",
    "import plotly.graph_objects as go\n",
    "from plotly.offline import init_notebook_mode\n",
    "init_notebook_mode() # Set notebook mode to work in offline. Might not be needed in your browser.\n",
    "\n",
    "x_0 = 0\n",
    "x_final = 10\n",
    "y_0 = 1\n",
    "h = 0.25\n",
    "\n",
    "def dx_dy(x, y):\n",
    "    return 2*math.cos(x)*y # right-hand side of the differential equation: 2(cos x)y\n",
    "\n",
    "n = (x_final - x_0)/h + 1\n",
    "print(f'There are n={int(n)} steps from {x_0} to {x_final} when h={h}.')\n",
    "\n",
    "# make sure we come up with an integer number for n\n",
    "assert n==int(n), 'hold on there, cowpoke! expecting an integer number of iterations for n'\n",
    "\n",
    "# Build arrays using the power of NumPy\n",
    "x = np.linspace(x_0, x_final, int(n))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4c569e3a",
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "# display values of array x (you can comment this out or delete it if you want)\n",
    "x"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c807e033",
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "y_euler = np.zeros(len(x))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f7cefd86",
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "# display initial values (all zeros!) of array y_euler (you can comment this out or delete it if you want)\n",
    "y_euler"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "dc6f7ec1",
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "y_euler[0] = y_0 # set the first (zeroeth) value of array y_euler to y_0\n",
    "\n",
    "# Euler's Method\n",
    "for i in range(1, len(x)): # iterate through x\n",
    "    y_euler[i] = y_euler[i-1] + h*dx_dy(x[i-1], y_euler[i-1]) # Euler's method is y_{n+1}=y_n+hf(x_n,y_n)\n",
    "#     print(f'iteration {i:2} x: {x[i]:.2f} y: {y_euler[i]:.4f}')\n",
    "# print()\n",
    "\n",
    "# print table\n",
    "assert len(x) == len(y_euler), \\\n",
    "    'whoopsie daisy! got different array lengths somewhere.'\n",
    "\n",
    "print('h =', h)\n",
    "print('n\\tx\\ty_euler') # \\t is a tab character\n",
    "\n",
    "for i in range(len(x)):\n",
    "    print(f\"{i}\\t{x[i]:.2f}\\t{y_euler[i]:.4f}\")\n",
    "    sleep(0.01) # pausing 1/100th of a second avoids a Jupyter Notebook error that sometimes happens when writing output too fast. If we printed from n=0 to n=500, it would only take 5 seconds. "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c8a058ea",
   "metadata": {},
   "source": [
    "### Note\n",
    "The following code is just a placeholder to demonstrate graphing two solution curves on the same axes. The placeholder code just adds 0.5 to each value of the Euler's method approximation. You should implement RK4. See  Zill p. 389, equation (6). "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a4f9871c",
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "# RK4 Method\n",
    "y_RK4 = np.zeros(len(x))\n",
    "y_RK4[0] = y_0  # Set the first (zeroeth) value of array y_RK4 to y_0\n",
    "\n",
    "# RK4 Method Calculation\n",
    "for i in range(1, len(x)):\n",
    "    k1 = h * dx_dy(x[i - 1], y_RK4[i - 1])\n",
    "    k2 = h * dx_dy(x[i - 1] + h / 2, y_RK4[i - 1] + k1 / 2)\n",
    "    k3 = h * dx_dy(x[i - 1] + h / 2, y_RK4[i - 1] + k2 / 2)\n",
    "    k4 = h * dx_dy(x[i - 1] + h, y_RK4[i - 1] + k3)\n",
    "    y_RK4[i] = y_RK4[i -\n",
    "                     1] + (k1 + 2 * k2 + 2 * k3 + k4) / 6 "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bae0e884",
   "metadata": {},
   "source": [
    "Now let's plot our approximate solutions. Put your name in place of \"Hi, my name is\" on the y-axis before turning in your plot."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a39af606",
   "metadata": {
    "scrolled": false,
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "# plot Euler's and RK4 approximations\n",
    "fig = go.Figure()\n",
    "fig.add_trace(go.Scatter(x=x, y=y_euler, mode='lines+markers', name=f\"Euler's method\"))\n",
    "fig.add_trace(go.Scatter(x=x, y=y_RK4, mode='lines+markers', name=f'4th-order Runge-Kutte (RK4)'))\n",
    "\n",
    "fig.update_layout(\n",
    "#     autosize=False,\n",
    "#     width=500,\n",
    "    height=800,\n",
    "    margin=dict(\n",
    "        l=50,\n",
    "        r=50,\n",
    "        b=100,\n",
    "        t=100,\n",
    "        pad=4\n",
    "    ),\n",
    "    title=rf\"$\\text{{Approximation of Initial-Value Problem }}y'(x)=2(\\cos x)y, y(0)=1\\text{{ using }}h={h}$\",\n",
    "    title_x=0.5,\n",
    "    xaxis_title=f'x-axis<br>FVS, Mr. Shadle, H Diff Eq, Oct 16, 2023',\n",
    "    yaxis_title=f'Hi, my name is<br>y-axis',\n",
    "    font=dict(\n",
    "        family=\"Courier New, monospace\",\n",
    "        size=18,\n",
    "        color=\"indigo\"\n",
    "    ),\n",
    "    paper_bgcolor=\"LightSteelBlue\",\n",
    "    legend=dict(\n",
    "        x=.35, #position of label on x-axis from 0 to 100%\n",
    "        y=1, #position of label on y-axis from 0 to 100%\n",
    "        traceorder=\"reversed\", # \"reversed\" puts RK4 first in the legend. set to \"normal\" to show Euler's first\n",
    "        font=dict(\n",
    "            family=\"sans-serif\",\n",
    "            size=12,\n",
    "            color=\"black\"\n",
    "        ),\n",
    "    )\n",
    ")\n",
    "fig.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "172d68a3",
   "metadata": {},
   "source": [
    "Now let's run it again for $h=0.1$ and $h=0.05$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cd8607bb",
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "R",
   "language": "R",
   "name": "ir"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "R",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.11.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}