nma2-metoda-strelby.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 76,
   "metadata": {},
   "outputs": [],
   "source": [
    "import math"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 77,
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_right_hand_side(x, y, dy) -> float: # returns right hand side of the initial problem\n",
    "    return 4*x*dy-(4*x*x-3)*y + pow(math.e, x*x)\n",
    "\n",
    "def get_coefficients(h, x, y, dy): # calculates the coefficients in Runge-Kutta-Nyström formula\n",
    "    f = get_right_hand_side(x, y, dy)\n",
    "    k1 = (h*h/2)*f\n",
    "    f = get_right_hand_side(x+h/2, y+h*dy/2+k1/4, dy+k1/h)\n",
    "    k2 = (h*h/2)*f\n",
    "    f = get_right_hand_side(x+h/2, y+h*dy/2+k2/4, dy+k2/h)\n",
    "    k3 = (h*h/2)*f\n",
    "    f = get_right_hand_side(x+h, y+h*dy+k3, dy+2*k3/h)\n",
    "    k4 = (h*h/2)*f    \n",
    "    return k1, k2, k3, k4\n",
    "\n",
    "def rkn_get_end(step, y0, dy0, start_time, stop_time): # returns the function value at the end of the interval\n",
    "    y = y0\n",
    "    dy = dy0\n",
    "    time = start_time\n",
    "    while(time <= stop_time):\n",
    "        k1, k2, k3, k4 = get_coefficients(step, time, y, dy)\n",
    "        y = y + step*dy + (1/3)*(k1+k2+k3)\n",
    "        dy = dy + (1/(3*step))*(k1+2*k2+2*k3+k4)\n",
    "        time += step\n",
    "    return y\n",
    "\n",
    "def get_init_dy(alpha, beta, step, start_time, stop_time): # returns initial values of the two derivations options\n",
    "    dy1 = -0.5\n",
    "    dy2 = 0.5\n",
    "    y1 = 0\n",
    "    y2 = 0\n",
    "    while((y1-beta)*(y2-beta) >= 0):\n",
    "        #print(\"dy1=\", dy1, \", dy2=\", dy2)\n",
    "        dy1 += dy1\n",
    "        dy2 += dy2\n",
    "        y1 = rkn_get_end(step, alpha, dy1, start_time, stop_time)\n",
    "        y2 = rkn_get_end(step, alpha, dy2, start_time, stop_time)\n",
    "        #print(\"y1=\", y1, \", y2=\", y2)\n",
    "    return dy1, dy2\n",
    "\n",
    "def rkn_shooting(epsilon, alpha, beta, step, start_time, stop_time): # returns dy sufficient for given epsilon\n",
    "    dy1, dy2 = get_init_dy(alpha, beta, step, start_time, stop_time)\n",
    "    y = beta + 2*epsilon # set far enough from the desired value\n",
    "    \n",
    "    while(abs(y-beta) > epsilon):\n",
    "        dy = (dy1+dy2)/2\n",
    "        #print(\"dy=(dy1+dy2)/2=(\",dy1,\"+\",dy2,\")/2=\",dy)\n",
    "        y = rkn_get_end(step, alpha, dy, start_time, stop_time)\n",
    "        if(y > beta):\n",
    "            dy2 = dy\n",
    "            #print(\"y=\", y, \", dy2=\", dy)\n",
    "        else:\n",
    "            dy1 = dy   \n",
    "            #print(\"y=\", y, \", dy1=\", dy)\n",
    "    return dy\n",
    "\n",
    "def rkn_solver(alpha, beta, step, epsilon, start_time, stop_time, output_path):\n",
    "    y = alpha\n",
    "    dy = rkn_shooting(epsilon, alpha, beta, step, start_time, stop_time)\n",
    "    time = start_time\n",
    "    f = open(output_path, \"w\") # overwrite content\n",
    "    f.write(\"#x\\t y\") # first line comment\n",
    "    f.write(\"\\n\")\n",
    "    f.close()\n",
    "    while(time <= stop_time):\n",
    "        f = open(output_path, \"a\") # now appending\n",
    "        f.write(str(round(time,3)))\n",
    "        f.write(\"\\t\")\n",
    "        f.write(str(round(y,3)))\n",
    "        f.write(\"\\n\")\n",
    "        f.close()\n",
    "        k1, k2, k3, k4 = get_coefficients(step, time, y, dy)\n",
    "        y = y + step*dy + (1/3)*(k1+k2+k3)\n",
    "        dy = dy + (1/(3*step))*(k1+2*k2+2*k3+k4)\n",
    "        time += step"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 94,
   "metadata": {},
   "outputs": [],
   "source": [
    "if __name__ == \"__main__\":\n",
    "    \n",
    "    # set global variables\n",
    "    alpha = 2 # value in 0\n",
    "    beta = 3 # value in stop_time\n",
    "    step = 0.01\n",
    "    start_time = 0\n",
    "    stop_time = 1 # corresponds with beta!\n",
    "    epsilon = 0.5\n",
    "    output_path = \"output.txt\"\n",
    "    \n",
    "    # shooting method\n",
    "    rkn_solver(alpha, beta, step, epsilon, start_time, stop_time, output_path)\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "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.7.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}