Partial Differential Equation.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sympy import *\n",
    "import numpy as np"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "\"\"\"\n",
    "author : YorozuyaSaint\n",
    "https://github.com/berylgithub\n",
    "\"\"\"\n",
    "\n",
    "def forward_diff_parabolic(u, range_x, range_t, h, k, alpha, log=True):\n",
    "    \"\"\"\n",
    "    numerical solution using forward difference method for parabolic heat equation\n",
    "    \"\"\"\n",
    "\n",
    "    x=Symbol('x')\n",
    "    t=Symbol('t')\n",
    "    if log:\n",
    "        print(\"u(x,t) = \",u(x,t))\n",
    "        print(\"h=\",h, \";x=\",range_x[1])\n",
    "        print(\"k=\",k, \";t=\",range_t[1])\n",
    "        print(\"alpha=\",alpha)\n",
    "        \n",
    "    l = (alpha**2)*k/(h**2)\n",
    "    max_iter_i = int((range_x[1]-range_x[0])/h) + 1\n",
    "    max_iter_j = int((range_t[1]-range_t[0])/k) + 1\n",
    "    m = max_iter_i-1\n",
    "    #initiate tridiagonal matrix A \n",
    "    A = np.zeros((m-1, m-1))\n",
    "    for i in range(m-1):\n",
    "        for j in range(m-1):\n",
    "            if i==0:\n",
    "                A[i] = np.array( [1-(2*l), l]+([0]*(len(A)-2)) )\n",
    "            elif i==len(A)-1:\n",
    "                A[i] = np.array( ([0]*(len(A)-2))+[l, 1-(2*l)] )\n",
    "            else:\n",
    "                A[i][i-1] = l\n",
    "                A[i][i] = (1-(2*l))\n",
    "                A[i][i+1] = l\n",
    "    #insert zeros for u(0,t) and u(l,t)\n",
    "    #A=np.insert(A, 0, np.zeros(m-1), 0)\n",
    "    #A=np.insert(A, len(A), np.zeros(m-1), 0)\n",
    "    if log:\n",
    "        print(\"A = \")\n",
    "        print(A)\n",
    "\n",
    "    #initiate vector x input\n",
    "    xi=np.array([h*i for i in range(max_iter_i)])\n",
    "\n",
    "    #initiate w with size (m+1) * (m+1) (2 rows for u(0,t) and u(l,t))\n",
    "    #w = np.zeros((max_iter_i, max_iter_j))\n",
    "    w = np.zeros((m-1, max_iter_j))\n",
    "    \n",
    "    #initiate vector w at j=0, where 1<=i<=m-1, where w(i,0) = f(xi) \n",
    "    w10 = np.array([u(x,t).evalf(subs={x:xi[i], t:0}) for i in range(1, m)])\n",
    "    \n",
    "    #add zeros at both 0 and l\n",
    "    #w10=np.insert(w10, 0, 0, 0)\n",
    "    #w10=np.insert(w10, len(w10), 0, 0)\n",
    "    \n",
    "    #set it as w(1), but transpose w first, then insert w(1,0)\n",
    "    w = np.transpose(w)\n",
    "    w[0] = w10\n",
    "    #w = np.transpose(w)\n",
    "    #print(np.transpose(w))\n",
    "    \n",
    "    #calculate the w at j>0 until j*k = t\n",
    "    for j in range(1, len(w)):\n",
    "        w[j] = np.matmul(A, np.transpose(w[j-1]))\n",
    "    if log :\n",
    "        print(\"w(\",len(w)-1,\") = \")\n",
    "        print(w[-1])\n",
    "    return w[-1]\n",
    "\n",
    "\n",
    "def backward_diff_parabolic(u, l, T, m, N, alpha, log=True):\n",
    "    \"\"\"\n",
    "    backward difference solution to parabolic heat partial differential equation\n",
    "    \"\"\"\n",
    "    x=Symbol('x')\n",
    "    t=Symbol('t')\n",
    "    if log:\n",
    "        print(\"u(x,t) = \",u(x,t))\n",
    "        print(\"l = \",l)\n",
    "        print(\"T = \",T)\n",
    "        print(\"m = \",m)\n",
    "        print(\"N = \",N)\n",
    "        print(\"alpha = \",alpha)\n",
    "        \n",
    "    #step 1\n",
    "    h=l/m\n",
    "    k=T/N\n",
    "    var_lambda = (alpha**2)*k/(h**2)\n",
    "    \n",
    "    #step 2\n",
    "    w = np.zeros(m)\n",
    "    for i in range(1, m):\n",
    "        w[i] = u(x,t).evalf(subs={x:(i*h), t:0})\n",
    "    \n",
    "    #step 3\n",
    "    l = np.zeros(m)\n",
    "    u = np.zeros(m)\n",
    "    l[1] = 1+(2*var_lambda)\n",
    "    u[1] = -var_lambda/l[1]\n",
    "    \n",
    "    #step 4\n",
    "    for i in range(2, m-1):\n",
    "        l[i] = 1 + (2*var_lambda) + (var_lambda*(u[i-1]))\n",
    "        u[i] = -var_lambda/l[i]\n",
    "        \n",
    "    #step 5\n",
    "    l[m-1] = 1 + (2*var_lambda) + (var_lambda*(u[m-2]))\n",
    "    \n",
    "    #step 6\n",
    "    if log :\n",
    "        print(\"Backward diff :\")\n",
    "    z = np.zeros(m)\n",
    "    for j in range(1, N+1):\n",
    "        #step 7\n",
    "        t = j*k\n",
    "        z[1] = w[1]/l[1]\n",
    "        #step 8\n",
    "        for i in range(2, m):\n",
    "            z[i] = (w[i] + var_lambda*z[i-1])/l[i]\n",
    "        #step 9\n",
    "        w[m-1] = z[m-1]\n",
    "        #step 10\n",
    "        for i in range(m-2, 0, -1):\n",
    "            w[i] = z[i] - (u[i]*w[i+1])\n",
    "        #step 11\n",
    "        if log:\n",
    "            print(\"j iter = \",j)\n",
    "            print(\"x\\tw[i]\")\n",
    "        for i in range(1, m):\n",
    "            x = i*h\n",
    "            \n",
    "            if log:\n",
    "                print(x,\"\\t\",w[i])\n",
    "        if log:\n",
    "            print()\n",
    "    \n",
    "    return w\n",
    "\n",
    "def crank_nicolson_parabolic(u, l, T, m, N, alpha, log=True):\n",
    "    \"\"\"\n",
    "    solves parabolic partial differential equation using crank-nicolson algorithm\n",
    "    \"\"\"\n",
    "    x=Symbol('x')\n",
    "    t=Symbol('t')\n",
    "    \n",
    "    if log:\n",
    "        print(\"u(x,t) = \",u(x,t))\n",
    "        print(\"l = \",l)\n",
    "        print(\"T = \",T)\n",
    "        print(\"m = \",m)\n",
    "        print(\"N = \",N)\n",
    "        print(\"alpha = \",alpha)\n",
    "        \n",
    "    #step 1\n",
    "    h=l/m\n",
    "    k=T/N\n",
    "    var_lambda = (alpha**2)*k/(h**2)\n",
    "    w = np.zeros(m+1)\n",
    "    \n",
    "    #step 2\n",
    "    for i in range(1, m):\n",
    "        w[i] = u(x,t).evalf(subs={x:(i*h), t:0})\n",
    "    \n",
    "    #step 3\n",
    "    l = np.zeros(m)\n",
    "    u = np.zeros(m)\n",
    "    l[1] = 1+var_lambda\n",
    "    u[1] = -var_lambda/(2*l[1])\n",
    "    \n",
    "    #step 4\n",
    "    for i in range(2, m-1):\n",
    "        l[i] = 1 + var_lambda + (var_lambda*(u[i-1])/2)\n",
    "        u[i] = -var_lambda/(2*l[i])\n",
    "        \n",
    "    #step 5\n",
    "    l[m-1] = 1 + var_lambda + (var_lambda*(u[m-2])/2)\n",
    "    \n",
    "    #step 6\n",
    "    if log :\n",
    "        print(\"Crank-Nicolson Algorithm :\")\n",
    "    z = np.zeros(m)\n",
    "    for j in range(1, N+1):\n",
    "        #step 7\n",
    "        t = j*k\n",
    "        z[1] = ( ((1-var_lambda)*w[1]) + (var_lambda*w[2]/2) )/l[1]\n",
    "        #step 8\n",
    "        for i in range(2, m):\n",
    "            z[i] = (((1-var_lambda)*w[i]) + (var_lambda*(w[i+1] + w[i-1] + z[i-1])/2) )/l[i]\n",
    "        #step 9\n",
    "        w[m-1] = z[m-1]\n",
    "        #step 10\n",
    "        for i in range(m-2, 0, -1):\n",
    "            w[i] = z[i] - (u[i]*w[i+1])\n",
    "        #step 11\n",
    "        if log:\n",
    "            print(\"j iter = \",j)\n",
    "            print(\"x\\tw[i]\")\n",
    "        for i in range(1, m):\n",
    "            x = i*h\n",
    "            \n",
    "            if log:\n",
    "                print(x,\"\\t\",w[i])\n",
    "        if log:\n",
    "            print()\n",
    "    \n",
    "    return w\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 235,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Example 1 page 727\n",
      "u(x,t) =  exp(-pi**2*t)*sin(pi*x)\n",
      "h= 0.1 ;x= 1\n",
      "k= 0.0005 ;t= 0.5\n",
      "alpha= 1\n",
      "A = \n",
      "[[0.9  0.05 0.   0.   0.   0.   0.   0.   0.  ]\n",
      " [0.05 0.9  0.05 0.   0.   0.   0.   0.   0.  ]\n",
      " [0.   0.05 0.9  0.05 0.   0.   0.   0.   0.  ]\n",
      " [0.   0.   0.05 0.9  0.05 0.   0.   0.   0.  ]\n",
      " [0.   0.   0.   0.05 0.9  0.05 0.   0.   0.  ]\n",
      " [0.   0.   0.   0.   0.05 0.9  0.05 0.   0.  ]\n",
      " [0.   0.   0.   0.   0.   0.05 0.9  0.05 0.  ]\n",
      " [0.   0.   0.   0.   0.   0.   0.05 0.9  0.05]\n",
      " [0.   0.   0.   0.   0.   0.   0.   0.05 0.9 ]]\n",
      "w( 1000 ) = \n",
      "[0.00228652 0.00434922 0.00598619 0.00703719 0.00739934 0.00703719\n",
      " 0.00598619 0.00434922 0.00228652]\n",
      "\n",
      "Result :\n",
      "x\tu(x,0.5)\tw(i,1000)\terror\n",
      "0.1 \t 0.00222241417851267 \t 0.0022865207865784528 \t 6.41066080657787e-5\n",
      "0.2 \t 0.00422728297276244 \t 0.004349220987439515 \t 0.000121938014677076\n",
      "0.30000000000000004 \t 0.00581835585642586 \t 0.005986189135245531 \t 0.000167833278819674\n",
      "0.4 \t 0.00683988752999332 \t 0.007037187382261514 \t 0.000197299852268192\n",
      "0.5 \t 0.00719188335582637 \t 0.007399336697334161 \t 0.000207453341507796\n",
      "0.6000000000000001 \t 0.00683988752999332 \t 0.007037187382261516 \t 0.000197299852268195\n",
      "0.7000000000000001 \t 0.00581835585642586 \t 0.005986189135245533 \t 0.000167833278819677\n",
      "0.8 \t 0.00422728297276244 \t 0.004349220987439515 \t 0.000121938014677077\n",
      "0.9 \t 0.00222241417851267 \t 0.0022865207865784528 \t 6.41066080657791e-5\n"
     ]
    }
   ],
   "source": [
    "#Example 1 page 727\n",
    "print(\"Example 1 page 727\")\n",
    "u = lambda x, t : exp((-1*(pi)**2)*t)*sin(pi*x)\n",
    "\n",
    "#variables initialization\n",
    "alpha = 1\n",
    "k = 0.0005\n",
    "h = 0.1\n",
    "rn_x = [0,1]\n",
    "rn_t = [0,0.5]\n",
    "x_iter = int((rn_x[1]-rn_x[0])/h)\n",
    "xi = np.array([h*i for i in range(1, x_iter)])\n",
    "\n",
    "#exact calculation of u(x,t)\n",
    "u_exact = np.array([u(x,t).evalf(subs={x:xi[i], t:rn_t[1]}) for i in range(len(xi))])\n",
    "\n",
    "#approximation of u(x,t)\n",
    "w_last = forward_diff_parabolic(u, rn_x, rn_t, h, k, alpha)\n",
    "\n",
    "print(\"\\nResult :\")\n",
    "print(\"x\\tu(x,0.5)\\tw(i,1000)\\terror\")\n",
    "for i in range(len(w_last)):\n",
    "    err = abs(u_exact[i] - w_last[i])\n",
    "    print(xi[i],\"\\t\",u_exact[i],\"\\t\",w_last[i],\"\\t\",err)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 44,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Example 1 page 731\n",
      "\n",
      "Result :\n",
      "x\tu(x,0.5)\tw(i,50)\terror\n",
      "0.314159265358979 \t 0.00600057911579518 \t 0.0028980166450535885 \t 0.00310256247074159\n",
      "0.628318530717959 \t 0.00661538012684316 \t 0.005512355229220069 \t 0.00110302489762309\n",
      "0.942477796076938 \t 0.00129259265614838 \t 0.0075871060767132344 \t 0.00629451342056486\n",
      "1.25663706143592 \t -0.00519035237459191 \t 0.008919178118941287 \t 0.0141095304935332\n",
      "1.57079632679490 \t -0.00701473268423438 \t 0.009378178863319288 \t 0.0163929115475537\n",
      "1.88495559215388 \t -0.00254308839819121 \t 0.008919178118941284 \t 0.0114622665171325\n",
      "2.19911485751286 \t 0.00421108721509700 \t 0.0075871060767132275 \t 0.00337601886161622\n",
      "2.51327412287183 \t 0.00718563074911635 \t 0.005512355229220063 \t 0.00167327551989628\n",
      "2.82743338823081 \t 0.00371076131639063 \t 0.0028980166450535854 \t 0.000812744671337047\n"
     ]
    }
   ],
   "source": [
    "#Example 2 page 731\n",
    "print(\"Example 1 page 731\")\n",
    "u = lambda x, t : exp((-1*(pi)**2)*t)*sin(pi*x)\n",
    "\n",
    "#variables initialization\n",
    "l = 1\n",
    "T = 0.5\n",
    "m = 10\n",
    "N = 50\n",
    "alpha = 1\n",
    "x_iter = int(m)\n",
    "xi = np.array([h*i for i in range(1, x_iter)])\n",
    "\n",
    "#exact calculation of u(x,t)\n",
    "u_exact = np.array([u(x,t).evalf(subs={x:xi[i], t:T}) for i in range(len(xi))])\n",
    "\n",
    "#approximation of u(x,t) using backward difference\n",
    "w_last = backward_diff_parabolic(u, l, T, m, N, alpha, log=False)\n",
    "w_last = w_last[1:]\n",
    "\n",
    "print(\"\\nResult :\")\n",
    "print(\"x\\tu(x,0.5)\\tw(i,50)\\terror\")\n",
    "for i in range(len(w_last)):\n",
    "    err = abs(u_exact[i] - w_last[i])\n",
    "    print(xi[i],\"\\t\",u_exact[i],\"\\t\",w_last[i],\"\\t\",err)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 46,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Example 3 page 735\n",
      "\n",
      "Result :\n",
      "x\tu(x,0.5)\tw(i,50)\terror\n",
      "0.314159265358979 \t 0.00600057911579518 \t 0.0023051233677887873 \t 0.00369545574800639\n",
      "0.628318530717959 \t 0.00661538012684316 \t 0.004384605199599512 \t 0.00223077492724365\n",
      "0.942477796076938 \t 0.00129259265614838 \t 0.0060348913251326705 \t 0.00474229866898429\n",
      "1.25663706143592 \t -0.00519035237459191 \t 0.007094440240201527 \t 0.0122847926147934\n",
      "1.57079632679490 \t -0.00701473268423438 \t 0.007459535914687764 \t 0.0144742685989221\n",
      "1.88495559215388 \t -0.00254308839819121 \t 0.007094440240201526 \t 0.00963752863839274\n",
      "2.19911485751286 \t 0.00421108721509700 \t 0.006034891325132668 \t 0.00182380411003566\n",
      "2.51327412287183 \t 0.00718563074911635 \t 0.00438460519959951 \t 0.00280102554951684\n",
      "2.82743338823081 \t 0.00371076131639063 \t 0.0023051233677887865 \t 0.00140563794860185\n"
     ]
    }
   ],
   "source": [
    "#Example 3 page 735\n",
    "print(\"Example 3 page 735\")\n",
    "u = lambda x, t : exp((-1*(pi)**2)*t)*sin(pi*x)\n",
    "\n",
    "#variables initialization\n",
    "l = 1\n",
    "T = 0.5\n",
    "m = 10\n",
    "N = 50\n",
    "alpha = 1\n",
    "x_iter = int(m)\n",
    "xi = np.array([h*i for i in range(1, x_iter)])\n",
    "\n",
    "#exact calculation of u(x,t)\n",
    "u_exact = np.array([u(x,t).evalf(subs={x:xi[i], t:T}) for i in range(len(xi))])\n",
    "\n",
    "#approximation of u(x,t) using backward difference\n",
    "w_last = crank_nicolson_parabolic(u, l, T, m, N, alpha, log=False)\n",
    "w_last = w_last[1:-1]\n",
    "\n",
    "print(\"\\nResult :\")\n",
    "print(\"x\\tu(x,0.5)\\tw(i,50)\\terror\")\n",
    "for i in range(len(w_last)):\n",
    "    err = abs(u_exact[i] - w_last[i])\n",
    "    print(xi[i],\"\\t\",u_exact[i],\"\\t\",w_last[i],\"\\t\",err)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exercise 12.2 No. 5,7,9"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## FTCS"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### No. 5.a"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Exercise 5.a - part 1\n",
      "u(x,t) =  exp(-4*pi**2*t)*sin(2*pi*x)\n",
      "h= 0.4 ;x= 2\n",
      "k= 0.1 ;t= 0.5\n",
      "alpha= 1\n",
      "A = \n",
      "[[-0.25   0.625  0.     0.   ]\n",
      " [ 0.625 -0.25   0.625  0.   ]\n",
      " [ 0.     0.625 -0.25   0.625]\n",
      " [ 0.     0.     0.625 -0.25 ]]\n",
      "w( 5 ) = \n",
      "[-1.87612234  3.03562971 -3.03562971  1.87612234]\n",
      "\n",
      "Exercise 5.a - part 2\n",
      "u(x,t) =  exp(-4*pi**2*t)*sin(2*pi*x)\n",
      "h= 0.4 ;x= 2\n",
      "k= 0.05 ;t= 0.5\n",
      "alpha= 1\n",
      "A = \n",
      "[[0.375  0.3125 0.     0.    ]\n",
      " [0.3125 0.375  0.3125 0.    ]\n",
      " [0.     0.3125 0.375  0.3125]\n",
      " [0.     0.     0.3125 0.375 ]]\n",
      "w( 10 ) = \n",
      "[ 8.50814426e-10 -1.37664664e-09  1.37664666e-09 -8.50814417e-10]\n",
      "[1.57249482678860e-9 -2.54435007687733e-9 2.54435007687733e-9\n",
      " -1.57249482678860e-9]\n",
      "\n",
      "Result :\n",
      "x\tu(x,0.5)\tw(i,5)\t|u(x,0.5)-w(i,5)|\tw(i,10)\t|u(x,0.5)-w(i,10)|\n",
      "0.4 \t 0.000000001572 \t -1.87612 \t 1.87612 \t 0.000000000851 \t 0.000000000722\n",
      "0.8 \t -0.000000002544 \t 3.03563 \t 3.03563 \t -0.000000001377 \t 0.000000001168\n",
      "1.2 \t 0.000000002544 \t -3.03563 \t 3.03563 \t 0.000000001377 \t 0.000000001168\n",
      "1.6 \t -0.000000001572 \t 1.87612 \t 1.87612 \t -0.000000000851 \t 0.000000000722\n"
     ]
    }
   ],
   "source": [
    "#Exercise 12.2, No.5,7,9\n",
    "u = lambda x,t: exp(-4*(pi**2)*t)*sin(2*pi*x)\n",
    "x = Symbol('x')\n",
    "t = Symbol('t')\n",
    "\n",
    "#init variable\n",
    "#No.5a part 1\n",
    "print(\"Exercise 5.a - part 1\")\n",
    "rn_x = [0,2]\n",
    "rn_t = [0,0.5]\n",
    "h = 0.4\n",
    "k = 0.1\n",
    "alpha = 1\n",
    "x_iter = int((rn_x[1]-rn_x[0])/h)\n",
    "xi = np.array([h*i for i in range(1, x_iter)])\n",
    "#exact calculation of u(x,t)\n",
    "u_exact = np.array([u(x,t).evalf(subs={x:xi[i], t:rn_t[1]}) for i in range(len(xi))])\n",
    "#approximation of u(x,t)\n",
    "w_last_1 = forward_diff_parabolic(u, rn_x, rn_t, h, k, alpha)\n",
    "\n",
    "#No.5a part 2\n",
    "print(\"\\nExercise 5.a - part 2\")\n",
    "h = 0.4\n",
    "k = 0.05\n",
    "x_iter = int((rn_x[1]-rn_x[0])/h)\n",
    "xi = np.array([h*i for i in range(1, x_iter)])\n",
    "#approximation of u(x,t)\n",
    "w_last_2 = forward_diff_parabolic(u, rn_x, rn_t, h, k, alpha)\n",
    "\n",
    "\n",
    "print(u_exact)\n",
    "\n",
    "print(\"\\nResult :\")\n",
    "print(\"x\\tu(x,0.5)\\tw(i,5)\\t|u(x,0.5)-w(i,5)|\\tw(i,10)\\t|u(x,0.5)-w(i,10)|\")\n",
    "for i in range(len(w_last)):\n",
    "    print(format(xi[i],'.1f'),\"\\t\",\n",
    "          format(u_exact[i], '.12f'),\"\\t\",\n",
    "          format(w_last_1[i], '.5f'),\"\\t\",\n",
    "          format(abs(u_exact[i] - w_last_1[i]), '.5f'),\"\\t\",\n",
    "          format(w_last_2[i], '.12f'),\"\\t\",\n",
    "          format(abs(u_exact[i] - w_last_2[i]), '.12f') )\n",
    "    \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### No. 5.b"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Exercise 5.b - part 1\n",
      "u(x,t) =  exp(-t)*sin(x)\n",
      "h= 0.314159265358979 ;x= 3.14159265358979\n",
      "k= 0.05 ;t= 0.5\n",
      "alpha= 1\n",
      "A = \n",
      "[[-0.01321184  0.50660592  0.          0.          0.          0.\n",
      "   0.          0.          0.        ]\n",
      " [ 0.50660592 -0.01321184  0.50660592  0.          0.          0.\n",
      "   0.          0.          0.        ]\n",
      " [ 0.          0.50660592 -0.01321184  0.50660592  0.          0.\n",
      "   0.          0.          0.        ]\n",
      " [ 0.          0.          0.50660592 -0.01321184  0.50660592  0.\n",
      "   0.          0.          0.        ]\n",
      " [ 0.          0.          0.          0.50660592 -0.01321184  0.50660592\n",
      "   0.          0.          0.        ]\n",
      " [ 0.          0.          0.          0.          0.50660592 -0.01321184\n",
      "   0.50660592  0.          0.        ]\n",
      " [ 0.          0.          0.          0.          0.          0.50660592\n",
      "  -0.01321184  0.50660592  0.        ]\n",
      " [ 0.          0.          0.          0.          0.          0.\n",
      "   0.50660592 -0.01321184  0.50660592]\n",
      " [ 0.          0.          0.          0.          0.          0.\n",
      "   0.          0.50660592 -0.01321184]]\n",
      "w( 10 ) = \n",
      "[0.18581972 0.35345011 0.48648234 0.57189429 0.60132525 0.57189429\n",
      " 0.48648234 0.35345011 0.18581972]\n",
      "[0.314159265358979 0.628318530717959 0.942477796076938 1.25663706143592\n",
      " 1.57079632679490 1.88495559215388 2.19911485751286 2.51327412287183\n",
      " 2.82743338823081]\n",
      "\n",
      "Result :\n",
      "x\tu(x,0.5)\tw(i,10)\t|u(x,0.5)-w(i,10)|\n",
      "0.31 \t 0.187428281461 \t 0.18582 \t 0.00161\n",
      "0.63 \t 0.356509776842 \t 0.35345 \t 0.00306\n",
      "0.94 \t 0.490693611317 \t 0.48648 \t 0.00421\n",
      "1.26 \t 0.576844936252 \t 0.57189 \t 0.00495\n",
      "1.57 \t 0.606530659713 \t 0.60133 \t 0.00521\n",
      "1.88 \t 0.576844936252 \t 0.57189 \t 0.00495\n",
      "2.20 \t 0.490693611317 \t 0.48648 \t 0.00421\n",
      "2.51 \t 0.356509776842 \t 0.35345 \t 0.00306\n",
      "2.83 \t 0.187428281461 \t 0.18582 \t 0.00161\n"
     ]
    }
   ],
   "source": [
    "#Exercise 12.2, No.5,7,9\n",
    "u = lambda x,t: exp(-t)*sin(x)\n",
    "x = Symbol('x')\n",
    "t = Symbol('t')\n",
    "\n",
    "#init variable\n",
    "#No.5a part 1\n",
    "print(\"Exercise 5.b - part 1\")\n",
    "rn_x = [0,pi.evalf()]\n",
    "rn_t = [0,0.5]\n",
    "h = pi.evalf()/10\n",
    "k = 0.05\n",
    "alpha = 1\n",
    "x_iter = int((rn_x[1]-rn_x[0])/h)\n",
    "xi = np.array([h*i for i in range(1, x_iter)])\n",
    "#exact calculation of u(x,t)\n",
    "u_exact = np.array([u(x,t).evalf(subs={x:xi[i], t:rn_t[1]}) for i in range(len(xi))])\n",
    "#approximation of u(x,t)\n",
    "w_last = forward_diff_parabolic(u, rn_x, rn_t, h, k, alpha)\n",
    "\n",
    "print(xi)\n",
    "print(\"\\nResult :\")\n",
    "print(\"x\\tu(x,0.5)\\tw(i,10)\\t|u(x,0.5)-w(i,10)|\")\n",
    "for i in range(len(w_last)):\n",
    "    print(format(xi[i],'.2f'),\"\\t\",\n",
    "          format(u_exact[i], '.12f'),\"\\t\",\n",
    "          format(w_last[i], '.5f'),\"\\t\",\n",
    "          format(abs(u_exact[i] - w_last[i]), '.5f') )\n",
    "    \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## BTCS"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### No 7.a"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Exercise 7.a - part 1\n",
      "u(x,t) =  exp(-4*pi**2*t)*sin(2*pi*x)\n",
      "l =  2\n",
      "T =  0.5\n",
      "m =  5\n",
      "N =  5\n",
      "alpha =  1\n",
      "Backward diff :\n",
      "j iter =  1\n",
      "x\tw[i]\n",
      "0.4 \t 0.18023194285349337\n",
      "0.8 \t -0.2916214093953809\n",
      "1.2000000000000002 \t 0.29162140939538106\n",
      "1.6 \t -0.18023194285349353\n",
      "\n",
      "j iter =  2\n",
      "x\tw[i]\n",
      "0.4 \t 0.055264321617551536\n",
      "0.8 \t -0.08941955074240389\n",
      "1.2000000000000002 \t 0.08941955074240394\n",
      "1.6 \t -0.055264321617551584\n",
      "\n",
      "j iter =  3\n",
      "x\tw[i]\n",
      "0.4 \t 0.016945637912425007\n",
      "0.8 \t -0.027418618103352436\n",
      "1.2000000000000002 \t 0.027418618103352464\n",
      "1.6 \t -0.01694563791242502\n",
      "\n",
      "j iter =  4\n",
      "x\tw[i]\n",
      "0.4 \t 0.00519602224100798\n",
      "0.8 \t -0.008407340592251284\n",
      "1.2000000000000002 \t 0.0084073405922513\n",
      "1.6 \t -0.005196022241007982\n",
      "\n",
      "j iter =  5\n",
      "x\tw[i]\n",
      "0.4 \t 0.0015932505620961882\n",
      "0.8 \t -0.0025779335620664907\n",
      "1.2000000000000002 \t 0.002577933562066501\n",
      "1.6 \t -0.001593250562096186\n",
      "\n",
      "\n",
      "Exercise 7.a - part 2\n",
      "u(x,t) =  exp(-4*pi**2*t)*sin(2*pi*x)\n",
      "l =  2\n",
      "T =  0.5\n",
      "m =  5\n",
      "N =  5\n",
      "alpha =  1\n",
      "Backward diff :\n",
      "j iter =  1\n",
      "x\tw[i]\n",
      "0.4 \t 0.18023194285349337\n",
      "0.8 \t -0.2916214093953809\n",
      "1.2000000000000002 \t 0.29162140939538106\n",
      "1.6 \t -0.18023194285349353\n",
      "\n",
      "j iter =  2\n",
      "x\tw[i]\n",
      "0.4 \t 0.055264321617551536\n",
      "0.8 \t -0.08941955074240389\n",
      "1.2000000000000002 \t 0.08941955074240394\n",
      "1.6 \t -0.055264321617551584\n",
      "\n",
      "j iter =  3\n",
      "x\tw[i]\n",
      "0.4 \t 0.016945637912425007\n",
      "0.8 \t -0.027418618103352436\n",
      "1.2000000000000002 \t 0.027418618103352464\n",
      "1.6 \t -0.01694563791242502\n",
      "\n",
      "j iter =  4\n",
      "x\tw[i]\n",
      "0.4 \t 0.00519602224100798\n",
      "0.8 \t -0.008407340592251284\n",
      "1.2000000000000002 \t 0.0084073405922513\n",
      "1.6 \t -0.005196022241007982\n",
      "\n",
      "j iter =  5\n",
      "x\tw[i]\n",
      "0.4 \t 0.0015932505620961882\n",
      "0.8 \t -0.0025779335620664907\n",
      "1.2000000000000002 \t 0.002577933562066501\n",
      "1.6 \t -0.001593250562096186\n",
      "\n",
      "\n",
      "Result :\n",
      "x\tu(x,0.5)\tw(i,5)\t|u(x,0.5)-w(i,5)|\tw(i,10)\t|u(x,0.5)-w(i,10)|\n",
      "0.4 \t 0.000000001572 \t 0.00159 \t 0.00159 \t 0.001593250562 \t 0.001593248990\n",
      "0.8 \t -0.000000002544 \t -0.00258 \t 0.00258 \t -0.002577933562 \t 0.002577931018\n",
      "1.2 \t 0.000000002544 \t 0.00258 \t 0.00258 \t 0.002577933562 \t 0.002577931018\n",
      "1.6 \t -0.000000001572 \t -0.00159 \t 0.00159 \t -0.001593250562 \t 0.001593248990\n"
     ]
    }
   ],
   "source": [
    "u = lambda x,t: exp(-4*(pi**2)*t)*sin(2*pi*x)\n",
    "x = Symbol('x')\n",
    "t = Symbol('t')\n",
    "\n",
    "#init variable\n",
    "#No.7a part 1\n",
    "print(\"Exercise 7.a - part 1\")\n",
    "l = 2\n",
    "T = 0.5\n",
    "h = 0.4\n",
    "k = 0.1\n",
    "m = int(l/h)\n",
    "N = int(T/k)\n",
    "alpha = 1\n",
    "x_iter = int(l/h)\n",
    "xi = np.array([h*i for i in range(1, x_iter)])\n",
    "#exact calculation of u(x,t)\n",
    "u_exact = np.array([u(x,t).evalf(subs={x:xi[i], t:T}) for i in range(len(xi))])\n",
    "\n",
    "#approximation of u(x,t) using backward difference\n",
    "w_last_1 = backward_diff_parabolic(u, l, T, m, N, alpha)\n",
    "w_last_1 = w_last_1[1:]\n",
    "\n",
    "#No.5a part 2\n",
    "print(\"\\nExercise 7.a - part 2\")\n",
    "h = 0.4\n",
    "k = 0.05\n",
    "x_iter = int(l/h)\n",
    "xi = np.array([h*i for i in range(1, x_iter)])\n",
    "#approximation of u(x,t)\n",
    "w_last_2 = backward_diff_parabolic(u, l, T, m, N, alpha)\n",
    "w_last_2 = w_last_2[1:]\n",
    "\n",
    "\n",
    "\n",
    "print(\"\\nResult :\")\n",
    "print(\"x\\tu(x,0.5)\\tw(i,5)\\t|u(x,0.5)-w(i,5)|\\tw(i,10)\\t|u(x,0.5)-w(i,10)|\")\n",
    "for i in range(len(w_last_1)):\n",
    "    print(format(xi[i],'.1f'),\"\\t\",\n",
    "          format(u_exact[i], '.12f'),\"\\t\",\n",
    "          format(w_last_1[i], '.5f'),\"\\t\",\n",
    "          format(abs(u_exact[i] - w_last_1[i]), '.5f'),\"\\t\",\n",
    "          format(w_last_2[i], '.12f'),\"\\t\",\n",
    "          format(abs(u_exact[i] - w_last_2[i]), '.12f') )\n",
    "    \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### No 7.b"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Exercise 7.b\n",
      "u(x,t) =  exp(-t)*sin(x)\n",
      "l =  3.14159265358979\n",
      "T =  0.5\n",
      "m =  10\n",
      "N =  10\n",
      "alpha =  1\n",
      "Backward diff :\n",
      "j iter =  1\n",
      "x\tw[i]\n",
      "0.314159265358979 \t 0.2944168293586503\n",
      "0.628318530717959 \t 0.5600140881370054\n",
      "0.942477796076938 \t 0.7707932661209245\n",
      "1.25663706143592 \t 0.9061218287844539\n",
      "1.57079632679490 \t 0.9527528735245483\n",
      "1.88495559215388 \t 0.906121828784454\n",
      "2.19911485751286 \t 0.7707932661209246\n",
      "2.51327412287183 \t 0.5600140881370055\n",
      "2.82743338823081 \t 0.29441682935865043\n",
      "\n",
      "j iter =  2\n",
      "x\tw[i]\n",
      "0.314159265358979 \t 0.28050648018544067\n",
      "0.628318530717959 \t 0.5335550316867614\n",
      "0.942477796076938 \t 0.7343754991900826\n",
      "1.25663706143592 \t 0.8633101761377072\n",
      "1.57079632679490 \t 0.9077380380092839\n",
      "1.88495559215388 \t 0.8633101761377073\n",
      "2.19911485751286 \t 0.7343754991900828\n",
      "2.51327412287183 \t 0.5335550316867617\n",
      "2.82743338823081 \t 0.2805064801854408\n",
      "\n",
      "j iter =  3\n",
      "x\tw[i]\n",
      "0.314159265358979 \t 0.2672533550389354\n",
      "0.628318530717959 \t 0.5083460896230434\n",
      "0.942477796076938 \t 0.6996783670993757\n",
      "1.25663706143592 \t 0.8225212510581844\n",
      "1.57079632679490 \t 0.864850024120881\n",
      "1.88495559215388 \t 0.8225212510581845\n",
      "2.19911485751286 \t 0.6996783670993759\n",
      "2.51327412287183 \t 0.5083460896230435\n",
      "2.82743338823081 \t 0.2672533550389355\n",
      "\n",
      "j iter =  4\n",
      "x\tw[i]\n",
      "0.314159265358979 \t 0.254626401972422\n",
      "0.628318530717959 \t 0.48432819763332213\n",
      "0.942477796076938 \t 0.666620574796894\n",
      "1.25663706143592 \t 0.7836594854806916\n",
      "1.57079632679490 \t 0.8239883456489442\n",
      "1.88495559215388 \t 0.7836594854806918\n",
      "2.19911485751286 \t 0.6666205747968942\n",
      "2.51327412287183 \t 0.48432819763332224\n",
      "2.82743338823081 \t 0.2546264019724221\n",
      "\n",
      "j iter =  5\n",
      "x\tw[i]\n",
      "0.314159265358979 \t 0.24259603615444175\n",
      "0.628318530717959 \t 0.46144508202411305\n",
      "0.942477796076938 \t 0.6351246681883268\n",
      "1.25663706143592 \t 0.7466338266564979\n",
      "1.57079632679490 \t 0.7850572640677702\n",
      "1.88495559215388 \t 0.7466338266564981\n",
      "2.19911485751286 \t 0.635124668188327\n",
      "2.51327412287183 \t 0.46144508202411305\n",
      "2.82743338823081 \t 0.24259603615444184\n",
      "\n",
      "j iter =  6\n",
      "x\tw[i]\n",
      "0.314159265358979 \t 0.23113407055180962\n",
      "0.628318530717959 \t 0.4396431278722446\n",
      "0.942477796076938 \t 0.6051168526627536\n",
      "1.25663706143592 \t 0.7113575238176079\n",
      "1.57079632679490 \t 0.7479655642218883\n",
      "1.88495559215388 \t 0.711357523817608\n",
      "2.19911485751286 \t 0.6051168526627538\n",
      "2.51327412287183 \t 0.4396431278722446\n",
      "2.82743338823081 \t 0.23113407055180962\n",
      "\n",
      "j iter =  7\n",
      "x\tw[i]\n",
      "0.314159265358979 \t 0.2202136498876623\n",
      "0.628318530717959 \t 0.4188712534056015\n",
      "0.942477796076938 \t 0.5765268201925693\n",
      "1.25663706143592 \t 0.6777479249205333\n",
      "1.57079632679490 \t 0.7126263406098141\n",
      "1.88495559215388 \t 0.6777479249205334\n",
      "2.19911485751286 \t 0.5765268201925694\n",
      "2.51327412287183 \t 0.4188712534056015\n",
      "2.82743338823081 \t 0.22021364988766232\n",
      "\n",
      "j iter =  8\n",
      "x\tw[i]\n",
      "0.314159265358979 \t 0.20980918771979906\n",
      "0.628318530717959 \t 0.399080790319016\n",
      "0.942477796076938 \t 0.549287584602441\n",
      "1.25663706143592 \t 0.6457262829933379\n",
      "1.57079632679490 \t 0.678956793765284\n",
      "1.88495559215388 \t 0.645726282993338\n",
      "2.19911485751286 \t 0.5492875846024411\n",
      "2.51327412287183 \t 0.39908079031901605\n",
      "2.82743338823081 \t 0.2098091877197991\n",
      "\n",
      "j iter =  9\n",
      "x\tw[i]\n",
      "0.314159265358979 \t 0.19989630649188994\n",
      "0.628318530717959 \t 0.38022536974489024\n",
      "0.942477796076938 \t 0.523335324621334\n",
      "1.25663706143592 \t 0.6152175716322283\n",
      "1.57079632679490 \t 0.6468780362588884\n",
      "1.88495559215388 \t 0.6152175716322285\n",
      "2.19911485751286 \t 0.5233353246213341\n",
      "2.51327412287183 \t 0.3802253697448903\n",
      "2.82743338823081 \t 0.19989630649188997\n",
      "\n",
      "j iter =  10\n",
      "x\tw[i]\n",
      "0.314159265358979 \t 0.19045178041709193\n",
      "0.628318530717959 \t 0.36226081361137796\n",
      "0.942477796076938 \t 0.49860923434987825\n",
      "1.25663706143592 \t 0.5861503092154002\n",
      "1.57079632679490 \t 0.6163149078655729\n",
      "1.88495559215388 \t 0.5861503092154003\n",
      "2.19911485751286 \t 0.49860923434987847\n",
      "2.51327412287183 \t 0.36226081361137813\n",
      "2.82743338823081 \t 0.190451780417092\n",
      "\n",
      "\n",
      "Result :\n",
      "x\tu(x,0.5)\tw(i,10)\t|u(x,0.5)-w(i,10)|\n",
      "0.31 \t 0.187428281461 \t 0.19045 \t 0.00302\n",
      "0.63 \t 0.356509776842 \t 0.36226 \t 0.00575\n",
      "0.94 \t 0.490693611317 \t 0.49861 \t 0.00792\n",
      "1.26 \t 0.576844936252 \t 0.58615 \t 0.00931\n",
      "1.57 \t 0.606530659713 \t 0.61631 \t 0.00978\n",
      "1.88 \t 0.576844936252 \t 0.58615 \t 0.00931\n",
      "2.20 \t 0.490693611317 \t 0.49861 \t 0.00792\n",
      "2.51 \t 0.356509776842 \t 0.36226 \t 0.00575\n",
      "2.83 \t 0.187428281461 \t 0.19045 \t 0.00302\n"
     ]
    }
   ],
   "source": [
    "u = lambda x,t: exp(-t)*sin(x)\n",
    "x = Symbol('x')\n",
    "t = Symbol('t')\n",
    "\n",
    "#init variable\n",
    "#No. 7.b\n",
    "print(\"Exercise 7.b\")\n",
    "l = pi.evalf()\n",
    "T = 0.5\n",
    "h = pi.evalf()/10\n",
    "k = 0.05\n",
    "m = int(l/h)\n",
    "N = int(T/k)\n",
    "alpha = 1\n",
    "x_iter = int(l/h)\n",
    "xi = np.array([h*i for i in range(1, x_iter)])\n",
    "#exact calculation of u(x,t)\n",
    "u_exact = np.array([u(x,t).evalf(subs={x:xi[i], t:T}) for i in range(len(xi))])\n",
    "#approximation of u(x,t)\n",
    "w_last = backward_diff_parabolic(u, l, T, m, N, alpha)\n",
    "w_last = w_last[1:]\n",
    "\n",
    "print(\"\\nResult :\")\n",
    "print(\"x\\tu(x,0.5)\\tw(i,10)\\t|u(x,0.5)-w(i,10)|\")\n",
    "for i in range(len(w_last)):\n",
    "    print(format(xi[i],'.2f'),\"\\t\",\n",
    "          format(u_exact[i], '.12f'),\"\\t\",\n",
    "          format(w_last[i], '.5f'),\"\\t\",\n",
    "          format(abs(u_exact[i] - w_last[i]), '.5f') )\n",
    "    \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Crank Nicolson"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### No 9.a"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Exercise 9.a - part 1\n",
      "u(x,t) =  exp(-4*pi**2*t)*sin(2*pi*x)\n",
      "l =  2\n",
      "T =  0.5\n",
      "m =  5\n",
      "N =  5\n",
      "alpha =  1\n",
      "Crank-Nicolson Algorithm :\n",
      "j iter =  1\n",
      "x\tw[i]\n",
      "0.4 \t -0.03603886603522751\n",
      "0.8 \t 0.05831211016100237\n",
      "1.2000000000000002 \t -0.058312110161002295\n",
      "1.6 \t 0.03603886603522746\n",
      "\n",
      "j iter =  2\n",
      "x\tw[i]\n",
      "0.4 \t 0.0022096503102783283\n",
      "0.8 \t -0.003575289305282028\n",
      "1.2000000000000002 \t 0.0035752893052820567\n",
      "1.6 \t -0.0022096503102783196\n",
      "\n",
      "j iter =  3\n",
      "x\tw[i]\n",
      "0.4 \t -0.0001354802476010191\n",
      "0.8 \t 0.0002192116454227337\n",
      "1.2000000000000002 \t -0.0002192116454227141\n",
      "1.6 \t 0.00013548024760103036\n",
      "\n",
      "j iter =  4\n",
      "x\tw[i]\n",
      "0.4 \t 8.306697853806942e-06\n",
      "0.8 \t -1.3440519461714596e-05\n",
      "1.2000000000000002 \t 1.344051946172977e-05\n",
      "1.6 \t -8.306697853797659e-06\n",
      "\n",
      "j iter =  5\n",
      "x\tw[i]\n",
      "0.4 \t -5.09308408092693e-07\n",
      "0.8 \t 8.240783150642584e-07\n",
      "1.2000000000000002 \t -8.240783150523501e-07\n",
      "1.6 \t 5.093084081000433e-07\n",
      "\n",
      "\n",
      "Exercise 9.a - part 2\n",
      "u(x,t) =  exp(-4*pi**2*t)*sin(2*pi*x)\n",
      "l =  2\n",
      "T =  0.5\n",
      "m =  5\n",
      "N =  5\n",
      "alpha =  1\n",
      "Crank-Nicolson Algorithm :\n",
      "j iter =  1\n",
      "x\tw[i]\n",
      "0.4 \t -0.03603886603522751\n",
      "0.8 \t 0.05831211016100237\n",
      "1.2000000000000002 \t -0.058312110161002295\n",
      "1.6 \t 0.03603886603522746\n",
      "\n",
      "j iter =  2\n",
      "x\tw[i]\n",
      "0.4 \t 0.0022096503102783283\n",
      "0.8 \t -0.003575289305282028\n",
      "1.2000000000000002 \t 0.0035752893052820567\n",
      "1.6 \t -0.0022096503102783196\n",
      "\n",
      "j iter =  3\n",
      "x\tw[i]\n",
      "0.4 \t -0.0001354802476010191\n",
      "0.8 \t 0.0002192116454227337\n",
      "1.2000000000000002 \t -0.0002192116454227141\n",
      "1.6 \t 0.00013548024760103036\n",
      "\n",
      "j iter =  4\n",
      "x\tw[i]\n",
      "0.4 \t 8.306697853806942e-06\n",
      "0.8 \t -1.3440519461714596e-05\n",
      "1.2000000000000002 \t 1.344051946172977e-05\n",
      "1.6 \t -8.306697853797659e-06\n",
      "\n",
      "j iter =  5\n",
      "x\tw[i]\n",
      "0.4 \t -5.09308408092693e-07\n",
      "0.8 \t 8.240783150642584e-07\n",
      "1.2000000000000002 \t -8.240783150523501e-07\n",
      "1.6 \t 5.093084081000433e-07\n",
      "\n",
      "\n",
      "Result :\n",
      "x\tu(x,0.5)\tw(i,5)\t|u(x,0.5)-w(i,5)|\tw(i,10)\t|u(x,0.5)-w(i,10)|\n",
      "0.4 \t 0.000000001572 \t -0.0000005 \t 0.0000005 \t -0.000000509308 \t 0.000000510881\n",
      "0.8 \t -0.000000002544 \t 0.0000008 \t 0.0000008 \t 0.000000824078 \t 0.000000826623\n",
      "1.2 \t 0.000000002544 \t -0.0000008 \t 0.0000008 \t -0.000000824078 \t 0.000000826623\n",
      "1.6 \t -0.000000001572 \t 0.0000005 \t 0.0000005 \t 0.000000509308 \t 0.000000510881\n"
     ]
    }
   ],
   "source": [
    "u = lambda x,t: exp(-4*(pi**2)*t)*sin(2*pi*x)\n",
    "x = Symbol('x')\n",
    "t = Symbol('t')\n",
    "\n",
    "#init variable\n",
    "#No.9a part 1\n",
    "print(\"Exercise 9.a - part 1\")\n",
    "l = 2\n",
    "T = 0.5\n",
    "h = 0.4\n",
    "k = 0.1\n",
    "m = int(l/h)\n",
    "N = int(T/k)\n",
    "alpha = 1\n",
    "x_iter = int(l/h)\n",
    "xi = np.array([h*i for i in range(1, x_iter)])\n",
    "#exact calculation of u(x,t)\n",
    "u_exact = np.array([u(x,t).evalf(subs={x:xi[i], t:T}) for i in range(len(xi))])\n",
    "\n",
    "#approximation of u(x,t) using backward difference\n",
    "w_last_1 = crank_nicolson_parabolic(u, l, T, m, N, alpha)\n",
    "w_last_1 = w_last_1[1:-1]\n",
    "\n",
    "#No.5a part 2\n",
    "print(\"\\nExercise 9.a - part 2\")\n",
    "h = 0.4\n",
    "k = 0.05\n",
    "x_iter = int(l/h)\n",
    "xi = np.array([h*i for i in range(1, x_iter)])\n",
    "#approximation of u(x,t)\n",
    "w_last_2 = crank_nicolson_parabolic(u, l, T, m, N, alpha)\n",
    "w_last_2 = w_last_2[1:-1]\n",
    "\n",
    "\n",
    "\n",
    "print(\"\\nResult :\")\n",
    "print(\"x\\tu(x,0.5)\\tw(i,5)\\t|u(x,0.5)-w(i,5)|\\tw(i,10)\\t|u(x,0.5)-w(i,10)|\")\n",
    "for i in range(len(w_last_1)):\n",
    "    print(format(xi[i],'.1f'),\"\\t\",\n",
    "          format(u_exact[i], '.12f'),\"\\t\",\n",
    "          format(w_last_1[i], '.7f'),\"\\t\",\n",
    "          format(abs(u_exact[i] - w_last_1[i]), '.7f'),\"\\t\",\n",
    "          format(w_last_2[i], '.12f'),\"\\t\",\n",
    "          format(abs(u_exact[i] - w_last_2[i]), '.12f') )\n",
    "    \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### No 9.b"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Exercise 9.b\n",
      "u(x,t) =  exp(-t)*sin(x)\n",
      "l =  3.14159265358979\n",
      "T =  0.5\n",
      "m =  10\n",
      "N =  10\n",
      "alpha =  1\n",
      "Crank-Nicolson Algorithm :\n",
      "j iter =  1\n",
      "x\tw[i]\n",
      "0.314159265358979 \t 0.294063576342132\n",
      "0.628318530717959 \t 0.559342160970484\n",
      "0.942477796076938 \t 0.7698684377170509\n",
      "1.25663706143592 \t 0.9050346277910579\n",
      "1.57079632679490 \t 0.951609722749838\n",
      "1.88495559215388 \t 0.9050346277910579\n",
      "2.19911485751286 \t 0.769868437717051\n",
      "2.51327412287183 \t 0.559342160970484\n",
      "2.82743338823081 \t 0.294063576342132\n",
      "\n",
      "j iter =  2\n",
      "x\tw[i]\n",
      "0.314159265358979 \t 0.27983375835376206\n",
      "0.628318530717959 \t 0.5322754387234175\n",
      "0.942477796076938 \t 0.7326142905697738\n",
      "1.25663706143592 \t 0.8612397512312513\n",
      "1.57079632679490 \t 0.9055610644320236\n",
      "1.88495559215388 \t 0.8612397512312515\n",
      "2.19911485751286 \t 0.7326142905697739\n",
      "2.51327412287183 \t 0.5322754387234175\n",
      "2.82743338823081 \t 0.2798337583537621\n",
      "\n",
      "j iter =  3\n",
      "x\tw[i]\n",
      "0.314159265358979 \t 0.2662925252030687\n",
      "0.628318530717959 \t 0.5065184826701397\n",
      "0.942477796076938 \t 0.6971628819316718\n",
      "1.25663706143592 \t 0.8195641208903106\n",
      "1.57079632679490 \t 0.861740713457206\n",
      "1.88495559215388 \t 0.8195641208903106\n",
      "2.19911485751286 \t 0.6971628819316718\n",
      "2.51327412287183 \t 0.5065184826701398\n",
      "2.82743338823081 \t 0.26629252520306873\n",
      "\n",
      "j iter =  4\n",
      "x\tw[i]\n",
      "0.314159265358979 \t 0.25340655607884643\n",
      "0.628318530717959 \t 0.4820079128614003\n",
      "0.942477796076938 \t 0.6634269767864762\n",
      "1.25663706143592 \t 0.7799051858561432\n",
      "1.57079632679490 \t 0.8200408414152596\n",
      "1.88495559215388 \t 0.7799051858561433\n",
      "2.19911485751286 \t 0.6634269767864762\n",
      "2.51327412287183 \t 0.4820079128614003\n",
      "2.82743338823081 \t 0.2534065560788465\n",
      "\n",
      "j iter =  5\n",
      "x\tw[i]\n",
      "0.314159265358979 \t 0.24114414257318234\n",
      "0.628318530717959 \t 0.4586834163212653\n",
      "0.942477796076938 \t 0.6313235614445418\n",
      "1.25663706143592 \t 0.7421653576837253\n",
      "1.57079632679490 \t 0.7803588377427189\n",
      "1.88495559215388 \t 0.7421653576837254\n",
      "2.19911485751286 \t 0.6313235614445418\n",
      "2.51327412287183 \t 0.4586834163212653\n",
      "2.82743338823081 \t 0.24114414257318237\n",
      "\n",
      "j iter =  6\n",
      "x\tw[i]\n",
      "0.314159265358979 \t 0.22947511065681347\n",
      "0.628318530717959 \t 0.43648759863542774\n",
      "0.942477796076938 \t 0.6007736392716807\n",
      "1.25663706143592 \t 0.7062517702599442\n",
      "1.57079632679490 \t 0.7425970572297346\n",
      "1.88495559215388 \t 0.7062517702599442\n",
      "2.19911485751286 \t 0.6007736392716808\n",
      "2.51327412287183 \t 0.43648759863542774\n",
      "2.82743338823081 \t 0.2294751106568135\n",
      "\n",
      "j iter =  7\n",
      "x\tw[i]\n",
      "0.314159265358979 \t 0.21837074643011867\n",
      "0.628318530717959 \t 0.41536584272120197\n",
      "0.942477796076938 \t 0.5717020363027354\n",
      "1.25663706143592 \t 0.6720760512886478\n",
      "1.57079632679490 \t 0.7066625797452333\n",
      "1.88495559215388 \t 0.6720760512886478\n",
      "2.19911485751286 \t 0.5717020363027354\n",
      "2.51327412287183 \t 0.415365842721202\n",
      "2.82743338823081 \t 0.2183707464301187\n",
      "\n",
      "j iter =  8\n",
      "x\tw[i]\n",
      "0.314159265358979 \t 0.2078037254670404\n",
      "0.628318530717959 \t 0.39526617443167583\n",
      "0.942477796076938 \t 0.5440372162615639\n",
      "1.25663706143592 \t 0.6395541048335962\n",
      "1.57079632679490 \t 0.6724669815890468\n",
      "1.88495559215388 \t 0.6395541048335961\n",
      "2.19911485751286 \t 0.5440372162615639\n",
      "2.51327412287183 \t 0.3952661744316758\n",
      "2.82743338823081 \t 0.2078037254670404\n",
      "\n",
      "j iter =  9\n",
      "x\tw[i]\n",
      "0.314159265358979 \t 0.19774804557807377\n",
      "0.628318530717959 \t 0.37613913466331617\n",
      "0.942477796076938 \t 0.5177111045322605\n",
      "1.25663706143592 \t 0.6086059043842194\n",
      "1.57079632679490 \t 0.6399261179083734\n",
      "1.88495559215388 \t 0.6086059043842192\n",
      "2.19911485751286 \t 0.5177111045322603\n",
      "2.51327412287183 \t 0.3761391346633161\n",
      "2.82743338823081 \t 0.19774804557807377\n",
      "\n",
      "j iter =  10\n",
      "x\tw[i]\n",
      "0.314159265358979 \t 0.18817896282687313\n",
      "0.628318530717959 \t 0.3579376576523224\n",
      "0.942477796076938 \t 0.49265892064845684\n",
      "1.25663706143592 \t 0.5791552959349814\n",
      "1.57079632679490 \t 0.6089599156431673\n",
      "1.88495559215388 \t 0.5791552959349813\n",
      "2.19911485751286 \t 0.4926589206484568\n",
      "2.51327412287183 \t 0.35793765765232227\n",
      "2.82743338823081 \t 0.1881789628268731\n",
      "\n",
      "\n",
      "Result :\n",
      "x\tu(x,0.5)\tw(i,10)\t|u(x,0.5)-w(i,10)|\n",
      "0.31 \t 0.187428281461 \t 0.18818 \t 0.00075\n",
      "0.63 \t 0.356509776842 \t 0.35794 \t 0.00143\n",
      "0.94 \t 0.490693611317 \t 0.49266 \t 0.00197\n",
      "1.26 \t 0.576844936252 \t 0.57916 \t 0.00231\n",
      "1.57 \t 0.606530659713 \t 0.60896 \t 0.00243\n",
      "1.88 \t 0.576844936252 \t 0.57916 \t 0.00231\n",
      "2.20 \t 0.490693611317 \t 0.49266 \t 0.00197\n",
      "2.51 \t 0.356509776842 \t 0.35794 \t 0.00143\n",
      "2.83 \t 0.187428281461 \t 0.18818 \t 0.00075\n"
     ]
    }
   ],
   "source": [
    "u = lambda x,t: exp(-t)*sin(x)\n",
    "x = Symbol('x')\n",
    "t = Symbol('t')\n",
    "\n",
    "#init variable\n",
    "#No. 7.b\n",
    "print(\"Exercise 9.b\")\n",
    "l = pi.evalf()\n",
    "T = 0.5\n",
    "h = pi.evalf()/10\n",
    "k = 0.05\n",
    "m = int(l/h)\n",
    "N = int(T/k)\n",
    "alpha = 1\n",
    "x_iter = int(l/h)\n",
    "xi = np.array([h*i for i in range(1, x_iter)])\n",
    "#exact calculation of u(x,t)\n",
    "u_exact = np.array([u(x,t).evalf(subs={x:xi[i], t:T}) for i in range(len(xi))])\n",
    "#approximation of u(x,t)\n",
    "w_last = crank_nicolson_parabolic(u, l, T, m, N, alpha)\n",
    "w_last = w_last[1:-1]\n",
    "\n",
    "print(\"\\nResult :\")\n",
    "print(\"x\\tu(x,0.5)\\tw(i,10)\\t|u(x,0.5)-w(i,10)|\")\n",
    "for i in range(len(w_last)):\n",
    "    print(format(xi[i],'.2f'),\"\\t\",\n",
    "          format(u_exact[i], '.12f'),\"\\t\",\n",
    "          format(w_last[i], '.5f'),\"\\t\",\n",
    "          format(abs(u_exact[i] - w_last[i]), '.5f') )"
   ]
  },
  {
   "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.6.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}