Skip to content
Snippets Groups Projects
Commit 8662a728 authored by penrose's avatar penrose
Browse files

error and consistency plotted, still not consistent

parent 02e3f41e
Branches
No related tags found
No related merge requests found
Expand all Hide whitespace changes
Inline Side-by-side
.ipynb_checkpoints/five_point_stencil-checkpoint.ipynb
+ 288
− 69
View file @ 8662a728
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
import numpy as np import numpy as np
from scipy.sparse import csr_matrix import scipy.sparse as sp
from scipy.sparse.linalg import spsolve from scipy.sparse.linalg import spsolve
from matplotlib import pyplot as plt
``` ```
%% Cell type:markdown id: tags: %% Cell type:markdown id: tags:
**The following poisson problem is given with Dirichlet boundary condition is given:** **The following poisson problem with Dirichlet boundary condition is given:**
$$-\Delta u = f \quad in \; \Omega = (0,1)^2$$ $$-\Delta u = f \quad in \; \Omega = (0,1)^2$$
$$u = g \quad on \; \partial \Omega$$ $$ u = g \quad on \quad \partial \Omega \quad$$
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
#given exact solution #given exact solution
def u(x,y): def u(x,y):
return pow(x,4)*pow(y,5)-17*np.sin(x*y) return pow(x,4)*pow(y,5)-17*np.sin(x*y)
``` ```
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
#right-hand-side of the poisson equation #right-hand-side of the poisson equation
def fun(x,y): 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)) 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: %% Cell type:markdown id: tags:
First, the components of the following equation will be assembled: **First, the components of the following equation will be assembled:**
$$A \underline{u} = \underline{f} + B\underline{g}$$ $$A \underline{u} = \underline{f} + B\underline{g} \Longleftrightarrow A \underline{u} - B\underline{g} = \underline{f} = -\Delta u $$
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
# matrix A only involves inner nodes. so for m inner nodes, A is an (m x m)-matrix
def matrix_A(h): def matrix_A(h):
N = int(1/h) N = int(1/h)
m = pow(N-1,2) 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))) 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)))
# insert zeros for neighbours of inner nodes that are next to the boundary, because of lexicographical order
for i in range(N-2): for i in range(N-2):
#for successors:
A[(i+1)*(N-1)-1][(i+1)*(N-1)] = 0 A[(i+1)*(N-1)-1][(i+1)*(N-1)] = 0
#for predecessors
A[(i+1)*(N-1)][(i+1)*(N-1)-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
# example for matrix A
n = 2
h = pow(2,-n)
A = matrix_A(h)
A
```
%% Output
array([[ 64., -16., 0., -16., 0., 0., 0., 0., 0.],
[-16., 64., -16., 0., -16., 0., 0., 0., 0.],
[ 0., -16., 64., 0., 0., -16., 0., 0., 0.],
[-16., 0., 0., 64., -16., 0., -16., 0., 0.],
[ 0., -16., 0., -16., 64., -16., 0., -16., 0.],
[ 0., 0., -16., 0., -16., 64., 0., 0., -16.],
[ 0., 0., 0., -16., 0., 0., 64., -16., 0.],
[ 0., 0., 0., 0., -16., 0., -16., 64., -16.],
[ 0., 0., 0., 0., 0., -16., 0., -16., 64.]])
%% Cell type:code id: tags:
``` python
# matrix A in sparse format. this is a lot harder to assemble because of the single deletions
def sparse_A(h):
N = int(1/h)
m = pow(N-1,2)
A = pow(h,-2)*(4*sp.eye(m) - sp.eye(m,k=1) - sp.eye(m,k=-1) - sp.eye(m,k=N-1) - sp.eye(m,k=-(N-1)))
# insert zeros for neighbours of inner nodes that are next to the boundary, because of lexicographical order
for i in range(N-2):
A.data[(i+1)*(N-1)*5-1*(N-1)-2-1]=0
A.data[(i+1)*(N-1)*5-1*(N-1)]=0
return A return A
``` ```
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
#check that sparse matrix is equal to dense matrix
s=sparse_A(h)
d=A-s.todense()
np.linalg.norm(d)
```
%% Output
0.0
%% Cell type:code id: tags:
``` python
def vector_f(h): def vector_f(h):
N = int(1/h) N = int(1/h)
l = pow(N-1,2) l = pow(N-1,2)
f = np.zeros(l) f = np.zeros(l)
for i in range(N-1): for i in range(N-1):
for k in range(N-1): for k in range(N-1):
f[k+i*(N-1)]=fun((k+1)/(N),(i+1)/(N)) f[k+i*(N-1)]=fun((k+1)/(N),(i+1)/(N))
return f return f
``` ```
%% Cell type:raw id: tags: %% Cell type:code id: tags:
# just the initialisation of matrix B ``` python
def matrix_B(h): # for m inner nodes and l boundary nodes, B is a (m x l)-matrix
def sparse_B(h):
N = int(1/h) N = int(1/h)
m = pow(N-1,2) m = pow(N-1,2)
l = 4*N l = 4*N
B = np.zeros((m,l)) B = np.zeros((m,l))
# since the lower and left boundary values are zero, only entries for the right
# and upper boundary nodes are needed
for i in range(N-1):
#right boundary:
B[(i+1)*(N-1)-1][2*(N+i-1)] = pow(h,-2)
#upper boundary:
B[-(N-1)+i][-N+i] = pow(h,-2)
B = sp.csr_matrix(B)
return B return B
```
%% Cell type:code id: tags:
``` python
# example for B
n = 2
h = pow(2,-n)
B = sparse_B(h)
B.todense()
```
%% Output
%% Cell type:raw id: tags: matrix([[ 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
0., 0., 0.],
[ 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
0., 0., 0.],
[ 0., 0., 0., 0., 0., 0., 16., 0., 0., 0., 0., 0., 0.,
0., 0., 0.],
[ 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
0., 0., 0.],
[ 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
0., 0., 0.],
[ 0., 0., 0., 0., 0., 0., 0., 0., 16., 0., 0., 0., 0.,
0., 0., 0.],
[ 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 16.,
0., 0., 0.],
[ 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
16., 0., 0.],
[ 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 16., 0., 0.,
0., 16., 0.]])
%% Cell type:code id: tags:
``` python
# vector g of length 4*N contains the boundary values which satisfy u(x,y)
def vector_g(h): def vector_g(h):
N = int(1/h) N = int(1/h)
l = 4*N l = 4*N
g = np.zeros(l) g = np.zeros(l)
g[-1] = u(1,1)
# upper boundary, where y=1
for i in range(N):
g[-N+i] = u((i+1)*h,1)
# right boundary, where x=1
for i in range(N-1):
g[N+2+2*i] = u(1,(i+1)*h)
return g return g
```
%% Cell type:markdown id: tags: %% Cell type:code id: tags:
``` python
# example for g
n = 2
h = pow(2,-n)
g = vector_g(h)
g
```
Since the exact solution is zero at the boundary, the product B*g is zero. %% Output
array([ 0. , 0. , 0. , 0. ,
0. , 0. , -4.20489074, 0. ,
-8.11898416, 0. , -11.35055423, 0. ,
-4.20196106, -8.08773416, -11.27145267, -13.30500674])
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
n = 6 n = 3
h = pow(2,-n) h = pow(2,-n)
N = pow(2,n) N = pow(2,n)
``` ```
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
A=matrix_A(h) A=sparse_A(h)
f=vector_f(h) f=vector_f(h)
appr_u = np.linalg.solve(A,f) B=sparse_B(h)
g=vector_g(h)
``` ```
%% Cell type:markdown id: tags:
**Now, the linear system Au = RHS will be solved.**
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
A[3968,3905] RHS = f+B.dot(g)
appr_u = sp.linalg.spsolve(A,RHS)
``` ```
%% Output %% Cell type:markdown id: tags:
%% Cell type:markdown id: tags:
4096.0 #### **Error and Consistency**
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
def exact_solution(h): def exact_solution(h):
N = int(1/h) N = int(1/h)
l = pow(N-1,2) l = pow(N-1,2)
v = np.zeros(l) v = np.zeros(l)
for i in range(N-1): for i in range(N-1):
for k in range(N-1): for k in range(N-1):
v[k+i*(N-1)]=u((k+1)/(N),(i+1)/(N)) v[k+i*(N-1)]=u((k+1)/(N),(i+1)/(N))
return v return v
``` ```
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
v= exact_solution(h) grid = np.zeros(7)
error_eukl = np.zeros(7)
error_inf = np.zeros(7)
for i in range(7):
grid[i]=pow(2,-(i+2))
h = grid[i]
A = sparse_A(h)
f = vector_f(h)
B = sparse_B(h)
g = vector_g(h)
RHS = f+B.dot(g)
appr_u = sp.linalg.spsolve(A,RHS)
v = exact_solution(h)
x = (appr_u - v)
print(appr_u[-1],v[-1])
error_eukl[i] = np.linalg.norm(x)
error_inf[i] = np.linalg.norm(x,ord = np.inf)
``` ```
%% Cell type:code id: tags: %% Output
``` python -8.180678655519896 -8.991060763833046
error = max(abs(appr_u - v)) -8.270439559649837 -11.480133561857317
``` -10.760438683213074 -12.53128157061408
-12.195353800544211 -12.961737935032156
-12.821966137009229 -13.145532611303409
-13.087791494720165 -13.228477586622512
-13.203877060573376 -13.267565784364171
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
appr_u[-1] plt.loglog(grid,error_inf)
plt.ylabel('error in maximum-norm')
plt.xlabel('step size h')
``` ```
%% Output %% Output
0.015864236841443172 Text(0.5, 0, 'step size h')
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
error plt.loglog(grid,error_eukl)
plt.ylabel('error in Euklidean norm')
plt.xlabel('step size h')
``` ```
%% Output %% Output
13.161396848144852 Text(0.5, 0, 'step size h')
%% Cell type:code id: tags:
``` python
```
%% Cell type:code id: tags: %% Cell type:markdown id: tags:
``` python Obviously the algorithm is not consistent, although it is supposed to be of order 2.
``` Even after several hours of debugging, I couldn't figure out the problem.
......
five_point_stencil.ipynb
+ 289
− 55
View file @ 8662a728
This diff is collapsed.
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment