diff --git a/five_point_stencil.ipynb b/five_point_stencil.ipynb
index 06119d92ee720e4ed76d3ce03e66ff8452a427c9..4dcdfe41985ef79b94637b08b9a7b9f8da96f610 100644
--- a/five_point_stencil.ipynb
+++ b/five_point_stencil.ipynb
@@ -1,10 +1,5 @@
{
"cells": [
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": []
- },
{
"cell_type": "code",
"execution_count": 1,
@@ -18,9 +13,9 @@
"cell_type": "markdown",
"metadata": {},
"source": [
- "First, the components of the following equation will be assembled:\n",
- "\n",
- "$$A \\underline{u} = \\underline{f} + B\\underline{g}$$"
+ "**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$$"
]
},
{
@@ -29,34 +24,34 @@
"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,
+ "execution_count": 85,
"metadata": {},
"outputs": [],
"source": [
- "def f(x,y):\n",
+ "#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))"
]
},
{
- "cell_type": "code",
- "execution_count": 18,
+ "cell_type": "markdown",
"metadata": {},
- "outputs": [],
"source": [
- "n = 3\n",
- "h = pow(2,-n)\n",
- "N = pow(2,n)"
+ "First, the components of the following equation will be assembled:\n",
+ "\n",
+ "$$A \\underline{u} = \\underline{f} + B\\underline{g}$$"
]
},
{
"cell_type": "code",
- "execution_count": 19,
+ "execution_count": 78,
"metadata": {},
"outputs": [],
"source": [
@@ -69,30 +64,27 @@
},
{
"cell_type": "code",
- "execution_count": 25,
+ "execution_count": 81,
"metadata": {},
- "outputs": [
- {
- "data": {
- "text/plain": [
- "(49, 49)"
- ]
- },
- "execution_count": 25,
- "metadata": {},
- "output_type": "execute_result"
- }
- ],
+ "outputs": [],
"source": [
- "matrix_A(h).shape"
+ "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": "code",
- "execution_count": 31,
+ "execution_count": 79,
"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",
@@ -103,7 +95,7 @@
},
{
"cell_type": "code",
- "execution_count": 1,
+ "execution_count": 80,
"metadata": {},
"outputs": [],
"source": [
@@ -115,66 +107,221 @@
" 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": 11,
+ "execution_count": 107,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "n = 3\n",
+ "h = pow(2,-n)\n",
+ "N = pow(2,n)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 108,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "A=matrix_A(h)\n",
+ "f=vector_f(h)\n",
+ "appr_u = np.linalg.solve(A,f)"
+ ]
+ },
+ {
+ "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,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
- "\u001b[0;31mSignature:\u001b[0m \u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0meye\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mN\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mM\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;32mNone\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mk\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mdtype\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;34m<\u001b[0m\u001b[0;32mclass\u001b[0m \u001b[0;34m'float'\u001b[0m\u001b[0;34m>\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0morder\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;34m'C'\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
- "\u001b[0;31mDocstring:\u001b[0m\n",
- "Return a 2-D array with ones on the diagonal and zeros elsewhere.\n",
- "\n",
- "Parameters\n",
- "----------\n",
- "N : int\n",
- " Number of rows in the output.\n",
- "M : int, optional\n",
- " Number of columns in the output. If None, defaults to `N`.\n",
- "k : int, optional\n",
- " Index of the diagonal: 0 (the default) refers to the main diagonal,\n",
- " a positive value refers to an upper diagonal, and a negative value\n",
- " to a lower diagonal.\n",
- "dtype : data-type, optional\n",
- " Data-type of the returned array.\n",
- "order : {'C', 'F'}, optional\n",
- " Whether the output should be stored in row-major (C-style) or\n",
- " column-major (Fortran-style) order in memory.\n",
- "\n",
- " .. versionadded:: 1.14.0\n",
- "\n",
- "Returns\n",
- "-------\n",
- "I : ndarray of shape (N,M)\n",
- " An array where all elements are equal to zero, except for the `k`-th\n",
- " diagonal, whose values are equal to one.\n",
- "\n",
- "See Also\n",
- "--------\n",
- "identity : (almost) equivalent function\n",
- "diag : diagonal 2-D array from a 1-D array specified by the user.\n",
- "\n",
- "Examples\n",
- "--------\n",
- ">>> np.eye(2, dtype=int)\n",
- "array([[1, 0],\n",
- " [0, 1]])\n",
- ">>> np.eye(3, k=1)\n",
- "array([[0., 1., 0.],\n",
- " [0., 0., 1.],\n",
- " [0., 0., 0.]])\n",
- "\u001b[0;31mFile:\u001b[0m ~/miniconda3/envs/Manhatten/lib/python3.7/site-packages/numpy/lib/twodim_base.py\n",
- "\u001b[0;31mType:\u001b[0m function\n"
+ "-13.228477586622512"
]
},
+ "execution_count": 91,
"metadata": {},
- "output_type": "display_data"
+ "output_type": "execute_result"
+ }
+ ],
+ "source": [
+ "u[-1]"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 92,
+ "metadata": {},
+ "outputs": [
+ {
+ "data": {
+ "text/plain": [
+ "0.025380345287435015"
+ ]
+ },
+ "execution_count": 92,
+ "metadata": {},
+ "output_type": "execute_result"
+ }
+ ],
+ "source": [
+ "appr_u[-1]"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 95,
+ "metadata": {},
+ "outputs": [
+ {
+ "data": {
+ "text/plain": [
+ "(16129,)"
+ ]
+ },
+ "execution_count": 95,
+ "metadata": {},
+ "output_type": "execute_result"
+ }
+ ],
+ "source": [
+ "f.shape"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 96,
+ "metadata": {},
+ "outputs": [
+ {
+ "data": {
+ "text/plain": [
+ "(16129, 16129)"
+ ]
+ },
+ "execution_count": 96,
+ "metadata": {},
+ "output_type": "execute_result"
+ }
+ ],
+ "source": [
+ "A.shape"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 97,
+ "metadata": {},
+ "outputs": [
+ {
+ "data": {
+ "text/plain": [
+ "(16129,)"
+ ]
+ },
+ "execution_count": 97,
+ "metadata": {},
+ "output_type": "execute_result"
+ }
+ ],
+ "source": [
+ "u.shape"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 99,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "ex=np.matmul(A,u)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 100,
+ "metadata": {},
+ "outputs": [
+ {
+ "data": {
+ "text/plain": [
+ "434807.0641872273"
+ ]
+ },
+ "execution_count": 100,
+ "metadata": {},
+ "output_type": "execute_result"
}
],
"source": [
- "np.eye?"
+ "ex[-1]"
]
},
{