all eq components assembled, still buggy

02e3f41e · penrose · 9c6c69bf · 02e3f41e
Commit 02e3f41e authored 4 years ago by penrose
--- a/five_point_stencil.ipynb
+++ b/five_point_stencil.ipynb
@@ -53,7 +53,7 @@
  },
  {
   "cell_type": "code",
-   "execution_count": 77,
+   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
@@ -70,7 +70,7 @@
  },
  {
   "cell_type": "code",
-   "execution_count": 78,
+   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
@@ -85,82 +85,67 @@
   ]
  },
  {
-   "cell_type": "raw",
+   "cell_type": "code",
+   "execution_count": 6,
   "metadata": {},
+   "outputs": [],
   "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",
+    "    for i in range(N-1):\n",
+    "        B[(i+1)*(N-1)-1][2*(N+i)]= pow(h,-2)\n",
+    "        B[(N-2)*(N-1)+i][N+2+2*i]=pow(h,-2)\n",
    "    return B"
   ]
  },
  {
-   "cell_type": "raw",
+   "cell_type": "code",
+   "execution_count": 7,
   "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",
+    "    for i in range(N):\n",
+    "        g[3*N+i] = u((i+1)*h,1)\n",
+    "    for i in range(N-1):    \n",
+    "        g[N+2+2*i] = u(1,(i+1)*h)\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": 115,
+   "execution_count": 64,
   "metadata": {},
   "outputs": [],
   "source": [
-    "n = 6\n",
+    "n = 2\n",
    "h = pow(2,-n)\n",
    "N = pow(2,n)"
   ]
  },
  {
   "cell_type": "code",
-   "execution_count": 116,
+   "execution_count": 65,
   "metadata": {},
   "outputs": [],
   "source": [
    "A=matrix_A(h)\n",
    "f=vector_f(h)\n",
-    "appr_u = np.linalg.solve(A,f)"
+    "B=matrix_B(h)\n",
-   ]
+    "g=vector_g(h)\n",
-  },
+    "RHS = f+np.matmul(B,g)\n",
-  {
+    "appr_u = np.linalg.solve(A,RHS)"
-   "cell_type": "code",
-   "execution_count": 129,
-   "metadata": {},
-   "outputs": [
-    {
-     "data": {
-      "text/plain": [
-       "4096.0"
-      ]
-     },
-     "execution_count": 129,
-     "metadata": {},
-     "output_type": "execute_result"
-    }
-   ],
-   "source": [
-    "A[3968,3905]"
   ]
  },
  {
   "cell_type": "code",
-   "execution_count": 130,
+   "execution_count": 66,
   "metadata": {},
   "outputs": [],
   "source": [
@@ -176,7 +161,7 @@
  },
  {
   "cell_type": "code",
-   "execution_count": 131,
+   "execution_count": 67,
   "metadata": {},
   "outputs": [],
   "source": [
@@ -185,7 +170,7 @@
  },
  {
   "cell_type": "code",
-   "execution_count": 132,
+   "execution_count": 68,
   "metadata": {},
   "outputs": [],
   "source": [
@@ -194,16 +179,16 @@
  },
  {
   "cell_type": "code",
-   "execution_count": 133,
+   "execution_count": 69,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
-       "0.015864236841443172"
+       "7.389346886114219"
      ]
     },
-     "execution_count": 133,
+     "execution_count": 69,
     "metadata": {},
     "output_type": "execute_result"
    }
@@ -214,16 +199,16 @@
  },
  {
   "cell_type": "code",
-   "execution_count": 134,
+   "execution_count": 70,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
-       "13.161396848144852"
+       "16.380407649947266"
      ]
     },
-     "execution_count": 134,
+     "execution_count": 70,
     "metadata": {},
     "output_type": "execute_result"
    }

 %% Cell type:code id: tags:
 ``` python
 import numpy as np
 from scipy.sparse import csr_matrix
 from scipy.sparse.linalg import spsolve
 ```
 %% Cell type:markdown id: tags:
 **The following poisson problem is given with Dirichlet boundary condition is given:**
 $$-\Delta u = f \quad in \; \Omega = (0,1)^2$$
 $$u = g \quad on \; \partial \Omega$$
 %% Cell type:code id: tags:
 ``` python
 #given exact solution
 def u(x,y):
    return pow(x,4)*pow(y,5)-17*np.sin(x*y)
 ```
 %% Cell type:code id: tags:
 ``` python
 #right-hand-side of the poisson equation
 def fun(x,y):
    return -(17*(x**2 + y**2)*np.sin(x*y) + 4*x**2*y**3+(5*x**2 + 3*y**2))
 ```
 %% Cell type:markdown id: tags:
 First, the components of the following equation will be assembled:
 $$A \underline{u} = \underline{f} + B\underline{g}$$
 %% Cell type:code id: tags:
 ``` python
 def matrix_A(h):
    N = int(1/h)
    m = pow(N-1,2)
    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)))
    for i in range(N-2):
        A[(i+1)*(N-1)-1][(i+1)*(N-1)] = 0
        A[(i+1)*(N-1)][(i+1)*(N-1)-1] = 0
    #A = csr_matrix(A)
    return A
 ```
 %% Cell type:code id: tags:
 ``` python
 def vector_f(h):
    N = int(1/h)
    l = pow(N-1,2)
    f = np.zeros(l)
    for i in range(N-1):
        for k in range(N-1):
            f[k+i*(N-1)]=fun((k+1)/(N),(i+1)/(N))
    return f
 ```
-%% Cell type:raw id: tags:
+%% Cell type:code id: tags:
-# just the initialisation of matrix B
+``` python
 def matrix_B(h):
    N = int(1/h)
    m = pow(N-1,2)
    l = 4*N
    B = np.zeros((m,l))
+    for i in range(N-1):
+        B[(i+1)*(N-1)-1][2*(N+i)]= pow(h,-2)
+        B[(N-2)*(N-1)+i][N+2+2*i]=pow(h,-2)
    return B
+```
-%% Cell type:raw id: tags:
+%% Cell type:code id: tags:
+``` python
 def vector_g(h):
    N = int(1/h)
    l = 4*N
    g = np.zeros(l)
-    g[-1] = u(1,1)
+    for i in range(N):
+        g[3*N+i] = u((i+1)*h,1)
+    for i in range(N-1):
+        g[N+2+2*i] = u(1,(i+1)*h)
    return g
+```
-%% Cell type:markdown id: tags:
-Since the exact solution is zero at the boundary, the product B*g is zero.
 %% Cell type:code id: tags:
 ``` python
-n = 6
+n = 2
 h = pow(2,-n)
 N = pow(2,n)
 ```
 %% Cell type:code id: tags:
 ``` python
 A=matrix_A(h)
 f=vector_f(h)
-appr_u = np.linalg.solve(A,f)
+B=matrix_B(h)
+g=vector_g(h)
+RHS = f+np.matmul(B,g)
+appr_u = np.linalg.solve(A,RHS)
 ```
 %% Cell type:code id: tags:
 ``` python
-A[3968,3905]
-```
-%% Output
-    4096.0
-%% Cell type:code id: tags:
-``` python
 def exact_solution(h):
    N = int(1/h)
    l = pow(N-1,2)
    v = np.zeros(l)
    for i in range(N-1):
        for k in range(N-1):
            v[k+i*(N-1)]=u((k+1)/(N),(i+1)/(N))
    return v
 ```
 %% Cell type:code id: tags:
 ``` python
 v= exact_solution(h)
 ```
 %% Cell type:code id: tags:
 ``` python
 error = max(abs(appr_u - v))
 ```
 %% Cell type:code id: tags:
 ``` python
 appr_u[-1]
 ```
 %% Output
-    0.015864236841443172
+    7.389346886114219
 %% Cell type:code id: tags:
 ``` python
 error
 ```
 %% Output
-    13.161396848144852
+    16.380407649947266
 %% Cell type:code id: tags:
 ``` python
 ```
 %% Cell type:code id: tags:
 ``` python
 ```