diff --git a/.ipynb_checkpoints/five_point_stencil-checkpoint.ipynb b/.ipynb_checkpoints/five_point_stencil-checkpoint.ipynb index 2846fb8afc67dae5e590104d3d7818bde7acedc7..099e014c0383328e9e37f5d7789f18358158f40b 100644 --- a/.ipynb_checkpoints/five_point_stencil-checkpoint.ipynb +++ b/.ipynb_checkpoints/five_point_stencil-checkpoint.ipynb @@ -1,17 +1,45 @@ { "cells": [ + { + "cell_type": "code", + "execution_count": 1, + "metadata": {}, + "outputs": [], + "source": [ + "import numpy as np\n", + "from scipy.sparse import csr_matrix\n", + "from scipy.sparse.linalg import spsolve" + ] + }, { "cell_type": "markdown", "metadata": {}, - "source": [] + "source": [ + "**The following poisson problem is given with Dirichlet boundary condition is given:**\n", + "$$-\\Delta u = f \\quad in \\; \\Omega = (0,1)^2$$\n", + "$$u = g \\quad on \\; \\partial \\Omega$$" + ] }, { "cell_type": "code", - "execution_count": 1, + "execution_count": 2, + "metadata": {}, + "outputs": [], + "source": [ + "#given exact solution\n", + "def u(x,y):\n", + " return pow(x,4)*pow(y,5)-17*np.sin(x*y)" + ] + }, + { + "cell_type": "code", + "execution_count": 3, "metadata": {}, "outputs": [], "source": [ - "import numpy as np" + "#right-hand-side of the poisson equation\n", + "def fun(x,y):\n", + " return -(17*(x**2 + y**2)*np.sin(x*y) + 4*x**2*y**3+(5*x**2 + 3*y**2))" ] }, { @@ -25,104 +53,192 @@ }, { "cell_type": "code", - "execution_count": 2, + "execution_count": 77, "metadata": {}, "outputs": [], "source": [ - "def u(x,y):\n", - " return pow(x,4)*pow(y,5)-17*np.sin(x*y)" + "def matrix_A(h):\n", + " N = int(1/h)\n", + " m = pow(N-1,2)\n", + " A = pow(h,-2)*(np.zeros((m,m))-4*np.eye(m)+np.eye(m,k=1)+np.eye(m,k=-1)+np.eye(m,k=N-1)+np.eye(m,k=-(N-1)))\n", + " for i in range(N-2):\n", + " A[(i+1)*(N-1)-1][(i+1)*(N-1)] = 0\n", + " A[(i+1)*(N-1)][(i+1)*(N-1)-1] = 0\n", + " #A = csr_matrix(A)\n", + " return A" ] }, { "cell_type": "code", - "execution_count": 13, + "execution_count": 78, "metadata": {}, "outputs": [], "source": [ - "def f(x,y):\n", - " return -(12*pow(x,2)*pow(y,5)+20*pow(x,4)*pow(y,3)+(pow(x,2)+pow(y,2))*17*np.sin(x*y))" + "def vector_f(h):\n", + " N = int(1/h)\n", + " l = pow(N-1,2)\n", + " f = np.zeros(l)\n", + " for i in range(N-1):\n", + " for k in range(N-1):\n", + " f[k+i*(N-1)]=fun((k+1)/(N),(i+1)/(N)) \n", + " return f" + ] + }, + { + "cell_type": "raw", + "metadata": {}, + "source": [ + "# just the initialisation of matrix B\n", + "def matrix_B(h):\n", + " N = int(1/h)\n", + " m = pow(N-1,2)\n", + " l = 4*N\n", + " B = np.zeros((m,l))\n", + " return B" + ] + }, + { + "cell_type": "raw", + "metadata": {}, + "source": [ + "def vector_g(h):\n", + " N = int(1/h)\n", + " l = 4*N\n", + " g = np.zeros(l)\n", + " g[-1] = u(1,1)\n", + " return g " + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Since the exact solution is zero at the boundary, the product B*g is zero." ] }, { "cell_type": "code", - "execution_count": 15, + "execution_count": 115, "metadata": {}, "outputs": [], "source": [ - "n = 2\n", + "n = 6\n", "h = pow(2,-n)\n", "N = pow(2,n)" ] }, { "cell_type": "code", - "execution_count": 29, + "execution_count": 116, "metadata": {}, "outputs": [], "source": [ - "def matrix_A(h):\n", - " N = int(1/h)\n", - " m = pow(N-1,2)\n", - " A = np.zeros((m,m))\n", - " return A" + "A=matrix_A(h)\n", + "f=vector_f(h)\n", + "appr_u = np.linalg.solve(A,f)" ] }, { "cell_type": "code", - "execution_count": 30, + "execution_count": 129, "metadata": {}, "outputs": [ { "data": { "text/plain": [ - "array([[0., 0., 0., 0., 0., 0., 0., 0., 0.],\n", - " [0., 0., 0., 0., 0., 0., 0., 0., 0.],\n", - " [0., 0., 0., 0., 0., 0., 0., 0., 0.],\n", - " [0., 0., 0., 0., 0., 0., 0., 0., 0.],\n", - " [0., 0., 0., 0., 0., 0., 0., 0., 0.],\n", - " [0., 0., 0., 0., 0., 0., 0., 0., 0.],\n", - " [0., 0., 0., 0., 0., 0., 0., 0., 0.],\n", - " [0., 0., 0., 0., 0., 0., 0., 0., 0.],\n", - " [0., 0., 0., 0., 0., 0., 0., 0., 0.]])" + "4096.0" ] }, - "execution_count": 30, + "execution_count": 129, "metadata": {}, "output_type": "execute_result" } ], "source": [ - "matrix_A(h)" + "A[3968,3905]" ] }, { "cell_type": "code", - "execution_count": 31, + "execution_count": 130, "metadata": {}, "outputs": [], "source": [ - "def matrix_B(h):\n", + "def exact_solution(h):\n", " N = int(1/h)\n", - " m = pow(N-1,2)\n", - " l = 4*N\n", - " B = np.zeros((m,l))\n", - " return B" + " l = pow(N-1,2)\n", + " v = np.zeros(l)\n", + " for i in range(N-1):\n", + " for k in range(N-1):\n", + " v[k+i*(N-1)]=u((k+1)/(N),(i+1)/(N)) \n", + " return v" ] }, { "cell_type": "code", - "execution_count": 1, + "execution_count": 131, "metadata": {}, "outputs": [], "source": [ - "def vector_g(h):\n", - " N = int(1/h)\n", - " l = 4*N\n", - " g = np.zeros(l)\n", - " g[-1] = u(1,1)\n", - " return g " + "v= exact_solution(h)" + ] + }, + { + "cell_type": "code", + "execution_count": 132, + "metadata": {}, + "outputs": [], + "source": [ + "error = max(abs(appr_u - v))" + ] + }, + { + "cell_type": "code", + "execution_count": 133, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "0.015864236841443172" + ] + }, + "execution_count": 133, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "appr_u[-1]" + ] + }, + { + "cell_type": "code", + "execution_count": 134, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "13.161396848144852" + ] + }, + "execution_count": 134, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "error\n" ] }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + }, { "cell_type": "code", "execution_count": null, diff --git a/five_point_stencil.ipynb b/five_point_stencil.ipynb index 4dcdfe41985ef79b94637b08b9a7b9f8da96f610..099e014c0383328e9e37f5d7789f18358158f40b 100644 --- a/five_point_stencil.ipynb +++ b/five_point_stencil.ipynb @@ -6,7 +6,9 @@ "metadata": {}, "outputs": [], "source": [ - "import numpy as np" + "import numpy as np\n", + "from scipy.sparse import csr_matrix\n", + "from scipy.sparse.linalg import spsolve" ] }, { @@ -31,13 +33,13 @@ }, { "cell_type": "code", - "execution_count": 85, + "execution_count": 3, "metadata": {}, "outputs": [], "source": [ "#right-hand-side of the poisson equation\n", "def fun(x,y):\n", - " return -(12*pow(x,2)*pow(y,5)+20*pow(x,4)*pow(y,3)+(pow(x,2)+pow(y,2))*17*np.sin(x*y))" + " return -(17*(x**2 + y**2)*np.sin(x*y) + 4*x**2*y**3+(5*x**2 + 3*y**2))" ] }, { @@ -51,7 +53,7 @@ }, { "cell_type": "code", - "execution_count": 78, + "execution_count": 77, "metadata": {}, "outputs": [], "source": [ @@ -59,12 +61,16 @@ " N = int(1/h)\n", " m = pow(N-1,2)\n", " A = pow(h,-2)*(np.zeros((m,m))-4*np.eye(m)+np.eye(m,k=1)+np.eye(m,k=-1)+np.eye(m,k=N-1)+np.eye(m,k=-(N-1)))\n", + " for i in range(N-2):\n", + " A[(i+1)*(N-1)-1][(i+1)*(N-1)] = 0\n", + " A[(i+1)*(N-1)][(i+1)*(N-1)-1] = 0\n", + " #A = csr_matrix(A)\n", " return A" ] }, { "cell_type": "code", - "execution_count": 81, + "execution_count": 78, "metadata": {}, "outputs": [], "source": [ @@ -79,10 +85,8 @@ ] }, { - "cell_type": "code", - "execution_count": 79, + "cell_type": "raw", "metadata": {}, - "outputs": [], "source": [ "# just the initialisation of matrix B\n", "def matrix_B(h):\n", @@ -94,10 +98,8 @@ ] }, { - "cell_type": "code", - "execution_count": 80, + "cell_type": "raw", "metadata": {}, - "outputs": [], "source": [ "def vector_g(h):\n", " N = int(1/h)\n", @@ -116,18 +118,18 @@ }, { "cell_type": "code", - "execution_count": 107, + "execution_count": 115, "metadata": {}, "outputs": [], "source": [ - "n = 3\n", + "n = 6\n", "h = pow(2,-n)\n", "N = pow(2,n)" ] }, { "cell_type": "code", - "execution_count": 108, + "execution_count": 116, "metadata": {}, "outputs": [], "source": [ @@ -138,192 +140,105 @@ }, { "cell_type": "code", - "execution_count": 109, - "metadata": {}, - "outputs": [], - "source": [ - "def exact_solution(h):\n", - " N = int(1/h)\n", - " l = pow(N-1,2)\n", - " v = np.zeros(l)\n", - " for i in range(N-1):\n", - " for k in range(N-1):\n", - " v[k+i*(N-1)]=u((k+1)/(N),(i+1)/(N)) \n", - " return v" - ] - }, - { - "cell_type": "code", - "execution_count": 110, - "metadata": {}, - "outputs": [ - { - "ename": "TypeError", - "evalue": "'numpy.ndarray' object is not callable", - "output_type": "error", - "traceback": [ - "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m", - "\u001b[0;31mTypeError\u001b[0m Traceback (most recent call last)", - "\u001b[0;32m<ipython-input-110-a5ff9eca95bf>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m\u001b[0m\n\u001b[0;32m----> 1\u001b[0;31m \u001b[0mu\u001b[0m\u001b[0;34m=\u001b[0m \u001b[0mexact_solution\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mh\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m", - "\u001b[0;32m<ipython-input-109-71aecd4cad10>\u001b[0m in \u001b[0;36mexact_solution\u001b[0;34m(h)\u001b[0m\n\u001b[1;32m 5\u001b[0m \u001b[0;32mfor\u001b[0m \u001b[0mi\u001b[0m \u001b[0;32min\u001b[0m \u001b[0mrange\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mN\u001b[0m\u001b[0;34m-\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 6\u001b[0m \u001b[0;32mfor\u001b[0m \u001b[0mk\u001b[0m \u001b[0;32min\u001b[0m \u001b[0mrange\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mN\u001b[0m\u001b[0;34m-\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m----> 7\u001b[0;31m \u001b[0mv\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mk\u001b[0m\u001b[0;34m+\u001b[0m\u001b[0mi\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mN\u001b[0m\u001b[0;34m-\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mu\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mk\u001b[0m\u001b[0;34m+\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m/\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mN\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mi\u001b[0m\u001b[0;34m+\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m/\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mN\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 8\u001b[0m \u001b[0;32mreturn\u001b[0m \u001b[0mv\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n", - "\u001b[0;31mTypeError\u001b[0m: 'numpy.ndarray' object is not callable" - ] - } - ], - "source": [ - "u= exact_solution(h)" - ] - }, - { - "cell_type": "code", - "execution_count": 106, - "metadata": {}, - "outputs": [ - { - "ename": "ValueError", - "evalue": "operands could not be broadcast together with shapes (49,) (16129,) ", - "output_type": "error", - "traceback": [ - "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m", - "\u001b[0;31mValueError\u001b[0m Traceback (most recent call last)", - "\u001b[0;32m<ipython-input-106-a6b4164c5b68>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m\u001b[0m\n\u001b[0;32m----> 1\u001b[0;31m \u001b[0mmax\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mabs\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mappr_u\u001b[0m \u001b[0;34m-\u001b[0m \u001b[0mu\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m", - "\u001b[0;31mValueError\u001b[0m: operands could not be broadcast together with shapes (49,) (16129,) " - ] - } - ], - "source": [ - "max(abs(appr_u - u))" - ] - }, - { - "cell_type": "code", - "execution_count": 91, + "execution_count": 129, "metadata": {}, "outputs": [ { "data": { "text/plain": [ - "-13.228477586622512" + "4096.0" ] }, - "execution_count": 91, + "execution_count": 129, "metadata": {}, "output_type": "execute_result" } ], "source": [ - "u[-1]" + "A[3968,3905]" ] }, { "cell_type": "code", - "execution_count": 92, + "execution_count": 130, "metadata": {}, - "outputs": [ - { - "data": { - "text/plain": [ - "0.025380345287435015" - ] - }, - "execution_count": 92, - "metadata": {}, - "output_type": "execute_result" - } - ], + "outputs": [], "source": [ - "appr_u[-1]" + "def exact_solution(h):\n", + " N = int(1/h)\n", + " l = pow(N-1,2)\n", + " v = np.zeros(l)\n", + " for i in range(N-1):\n", + " for k in range(N-1):\n", + " v[k+i*(N-1)]=u((k+1)/(N),(i+1)/(N)) \n", + " return v" ] }, { "cell_type": "code", - "execution_count": 95, + "execution_count": 131, "metadata": {}, - "outputs": [ - { - "data": { - "text/plain": [ - "(16129,)" - ] - }, - "execution_count": 95, - "metadata": {}, - "output_type": "execute_result" - } - ], + "outputs": [], "source": [ - "f.shape" + "v= exact_solution(h)" ] }, { "cell_type": "code", - "execution_count": 96, + "execution_count": 132, "metadata": {}, - "outputs": [ - { - "data": { - "text/plain": [ - "(16129, 16129)" - ] - }, - "execution_count": 96, - "metadata": {}, - "output_type": "execute_result" - } - ], + "outputs": [], "source": [ - "A.shape" + "error = max(abs(appr_u - v))" ] }, { "cell_type": "code", - "execution_count": 97, + "execution_count": 133, "metadata": {}, "outputs": [ { "data": { "text/plain": [ - "(16129,)" + "0.015864236841443172" ] }, - "execution_count": 97, + "execution_count": 133, "metadata": {}, "output_type": "execute_result" } ], "source": [ - "u.shape" - ] - }, - { - "cell_type": "code", - "execution_count": 99, - "metadata": {}, - "outputs": [], - "source": [ - "ex=np.matmul(A,u)" + "appr_u[-1]" ] }, { "cell_type": "code", - "execution_count": 100, + "execution_count": 134, "metadata": {}, "outputs": [ { "data": { "text/plain": [ - "434807.0641872273" + "13.161396848144852" ] }, - "execution_count": 100, + "execution_count": 134, "metadata": {}, "output_type": "execute_result" } ], "source": [ - "ex[-1]" + "error\n" ] }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + }, { "cell_type": "code", "execution_count": null,