Skip to content
Snippets Groups Projects

Compare revisions

Changes are shown as if the source revision was being merged into the target revision. Learn more about comparing revisions.

Source

Select target project
No results found
Select Git revision
  • master
1 result

Target

Select target project
  • felixhenneke/exercise_problems_03
  • numerics_3/exercise_problems_03
2 results
Select Git revision
  • master
1 result
Show changes
Commits on Source (3)
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
import numpy as np import numpy as np
import scipy.sparse as sp import scipy.sparse as sp
from scipy.sparse.linalg import spsolve from scipy.sparse.linalg import spsolve
from matplotlib import pyplot as plt from matplotlib import pyplot as plt
``` ```
%% Cell type:markdown id: tags: %% Cell type:markdown id: tags:
**The following poisson problem 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 \quad \partial \Omega \quad$$ $$ 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} \Longleftrightarrow A \underline{u} - B\underline{g} = \underline{f} = -\Delta u $$ $$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 # 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 # 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: #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 #for predecessors
A[(i+1)*(N-1)][(i+1)*(N-1)-1] = 0 A[(i+1)*(N-1)][(i+1)*(N-1)-1] = 0
return A return A
``` ```
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
# example for matrix A # example for matrix A
n = 2 n = 2
h = pow(2,-n) h = pow(2,-n)
A = matrix_A(h) A = matrix_A(h)
A A
``` ```
%% Output %% Output
array([[ 64., -16., 0., -16., 0., 0., 0., 0., 0.], array([[ 64., -16., 0., -16., 0., 0., 0., 0., 0.],
[-16., 64., -16., 0., -16., 0., 0., 0., 0.], [-16., 64., -16., 0., -16., 0., 0., 0., 0.],
[ 0., -16., 64., 0., 0., -16., 0., 0., 0.], [ 0., -16., 64., 0., 0., -16., 0., 0., 0.],
[-16., 0., 0., 64., -16., 0., -16., 0., 0.], [-16., 0., 0., 64., -16., 0., -16., 0., 0.],
[ 0., -16., 0., -16., 64., -16., 0., -16., 0.], [ 0., -16., 0., -16., 64., -16., 0., -16., 0.],
[ 0., 0., -16., 0., -16., 64., 0., 0., -16.], [ 0., 0., -16., 0., -16., 64., 0., 0., -16.],
[ 0., 0., 0., -16., 0., 0., 64., -16., 0.], [ 0., 0., 0., -16., 0., 0., 64., -16., 0.],
[ 0., 0., 0., 0., -16., 0., -16., 64., -16.], [ 0., 0., 0., 0., -16., 0., -16., 64., -16.],
[ 0., 0., 0., 0., 0., -16., 0., -16., 64.]]) [ 0., 0., 0., 0., 0., -16., 0., -16., 64.]])
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
# matrix A in sparse format. this is a lot harder to assemble because of the single deletions # matrix A in sparse format. this is a lot harder to assemble because of the single deletions
def sparse_A(h): def sparse_A(h):
N = int(1/h) N = int(1/h)
m = pow(N-1,2) 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))) 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 # 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):
A.data[(i+1)*(N-1)*5-1*(N-1)-2-1]=0 A.data[(i+1)*(N-1)*5-1*(N-1)-2-1]=0
A.data[(i+1)*(N-1)*5-1*(N-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 #check that sparse matrix is equal to dense matrix
s=sparse_A(h) s=sparse_A(h)
d=A-s.todense() d=A-s.todense()
np.linalg.norm(d) np.linalg.norm(d)
``` ```
%% Output %% Output
0.0 0.0
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` 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:code id: tags: %% Cell type:code id: tags:
``` python ``` python
# for m inner nodes and l boundary nodes, B is a (m x l)-matrix # for m inner nodes and l boundary nodes, B is a (m x l)-matrix
def sparse_B(h): 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 # since the lower and left boundary values are zero, only entries for the right
# and upper boundary nodes are needed # and upper boundary nodes are needed
for i in range(N-1): for i in range(N-1):
#right boundary: #right boundary:
B[(i+1)*(N-1)-1][2*(N+i-1)] = pow(h,-2) B[(i+1)*(N-1)-1][N + 2 * (i + 1)] = pow(h,-2)
#upper boundary: #upper boundary:
B[-(N-1)+i][-N+i] = pow(h,-2) B[-(N-1)+i][-N+i] = pow(h,-2)
B = sp.csr_matrix(B) B = sp.csr_matrix(B)
return B return B
``` ```
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
# example for B # example for B
n = 2 n = 2
h = pow(2,-n) h = pow(2,-n)
B = sparse_B(h) B = sparse_B(h)
B.todense() B.todense()
``` ```
%% Output %% Output
matrix([[ 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 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., 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., 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., 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., 16., 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., 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.,
16., 0., 0.], 16., 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.]]) 0., 16., 0.]])
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
# vector g of length 4*N contains the boundary values which satisfy u(x,y) # 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)
# upper boundary, where y=1 # upper boundary, where y=1
for i in range(N): for i in range(N):
g[-N+i] = u((i+1)*h,1) g[-N+i] = u((i+1)*h,1)
# right boundary, where x=1 # right boundary, where x=1
for i in range(N-1): for i in range(N-1):
g[N+2+2*i] = u(1,(i+1)*h) g[N+2+2*i] = u(1,(i+1)*h)
return g return g
``` ```
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
# example for g # example for g
n = 2 n = 2
h = pow(2,-n) h = pow(2,-n)
g = vector_g(h) g = vector_g(h)
g g
``` ```
%% Output %% Output
array([ 0. , 0. , 0. , 0. , array([ 0. , 0. , 0. , 0. ,
0. , 0. , -4.20489074, 0. , 0. , 0. , -4.20489074, 0. ,
-8.11898416, 0. , -11.35055423, 0. , -8.11898416, 0. , -11.35055423, 0. ,
-4.20196106, -8.08773416, -11.27145267, -13.30500674]) -4.20196106, -8.08773416, -11.27145267, -13.30500674])
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
n = 3 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=sparse_A(h) A=sparse_A(h)
f=vector_f(h) f=vector_f(h)
B=sparse_B(h) B=sparse_B(h)
g=vector_g(h) g=vector_g(h)
``` ```
%% Cell type:markdown id: tags: %% Cell type:markdown id: tags:
**Now, the linear system Au = RHS will be solved.** **Now, the linear system Au = RHS will be solved.**
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
RHS = f+B.dot(g) RHS = f+B.dot(g)
appr_u = sp.linalg.spsolve(A,RHS) appr_u = sp.linalg.spsolve(A,RHS)
``` ```
%% Cell type:markdown id: tags: %% Cell type:markdown id: tags:
%% Cell type:markdown id: tags:
#### **Error and Consistency** #### **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
grid = np.zeros(7) grid = np.zeros(7)
error_eukl = np.zeros(7) error_eukl = np.zeros(7)
error_inf = np.zeros(7) error_inf = np.zeros(7)
for i in range(7): for i in range(7):
grid[i]=pow(2,-(i+2)) grid[i]=pow(2,-(i+2))
h = grid[i] h = grid[i]
A = sparse_A(h) A = sparse_A(h)
f = vector_f(h) f = vector_f(h)
B = sparse_B(h) B = sparse_B(h)
g = vector_g(h) g = vector_g(h)
RHS = f+B.dot(g) RHS = f+B.dot(g)
appr_u = sp.linalg.spsolve(A,RHS) appr_u = sp.linalg.spsolve(A,RHS)
v = exact_solution(h) v = exact_solution(h)
x = (appr_u - v) x = (appr_u - v)
print(appr_u[-1],v[-1]) print(appr_u[-1],v[-1])
error_eukl[i] = np.linalg.norm(x) error_eukl[i] = np.linalg.norm(x)
error_inf[i] = np.linalg.norm(x,ord = np.inf) error_inf[i] = np.linalg.norm(x,ord = np.inf)
``` ```
%% Output %% Output
-8.180678655519896 -8.991060763833046 -8.987207596607227 -8.991060763833046
-8.270439559649837 -11.480133561857317 -11.479386388270894 -11.480133561857317
-10.760438683213074 -12.53128157061408 -12.531189312354407 -12.53128157061408
-12.195353800544211 -12.961737935032156 -12.96172877159427 -12.961737935032156
-12.821966137009229 -13.145532611303409 -13.145531806826055 -13.145532611303409
-13.087791494720165 -13.228477586622512 -13.228477521210696 -13.228477586622512
-13.203877060573376 -13.267565784364171 -13.267565779309312 -13.267565784364171
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
plt.loglog(grid,error_inf) plt.loglog(grid,error_inf)
plt.ylabel('error in maximum-norm') plt.ylabel('error in maximum-norm')
plt.xlabel('step size h') plt.xlabel('step size h')
``` ```
%% Output %% Output
Text(0.5, 0, 'step size h') Text(0.5, 0, 'step size h')
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
plt.loglog(grid,error_eukl) plt.loglog(grid,error_eukl)
plt.ylabel('error in Euklidean norm') plt.ylabel('error in Euklidean norm')
plt.xlabel('step size h') plt.xlabel('step size h')
``` ```
%% Output %% Output
Text(0.5, 0, 'step size h') Text(0.5, 0, 'step size h')
%% Cell type:markdown id: tags:
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.
......
%% Cell type:code id: tags:
``` python
import numpy as np
import scipy.sparse as sp
from scipy.sparse.linalg import spsolve
from matplotlib import pyplot as plt
```
%% Cell type:markdown id: tags:
**The following poisson problem with Dirichlet boundary condition is given:**
$$-\Delta u = f \quad in \; \Omega = (0,1)^2$$
$$ u = g \quad on \quad \partial \Omega \quad$$
%% 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} \Longleftrightarrow A \underline{u} - B\underline{g} = \underline{f} = -\Delta u $$
%% Cell type:code id: tags:
``` python
# matrix A only involves inner nodes. so for m inner nodes, A is an (m x m)-matrix
def sparse_A(h):
N = int(1/h)
m = pow(N-1,2)
A = np.zeros((m,m))+4*np.eye(m) - 2*np.eye(m,k=1) - 2*np.eye(m,k=-1) - 2*np.eye(m,k=N-1) - 2*np.eye(m,k=-(N-1))
A = A + np.eye(m,k=-N) + np.eye(m,k=-(N-2)) + np.eye(m,k=N-2) + np.eye(m,k=N)
A = pow(h,-4)*A
# insert zeros for neighbours of inner nodes that are next to the boundary, because of lexicographical order
for i in range(N-2):
#right boundary, successors:
A[(i+1)*(N-1)-1][(i+1)*(N-1)] = 0
#left boundary, predecessors:
A[(i+1)*(N-1)][(i+1)*(N-1)-1] = 0
#left boundary, diagonally below:
A[(i+1)*(N-1)][(i+1)*(N-1)-(N-2)] = 0
for i in range(N-1):
#right boundary:
A[(i+1)*(N-1)-1][(i+1)*(N-1)-1-2] = 0
#left boundary
A[i*(N-1)][i*(N-1)+(2)] = 0
#right boundary, diagonally below:
A[i*(N-1)+N-2][i*(N-1)] = 0
#left boundary, diagonally above:
A[i*(N-1)][i*(N-1)+N-2] = 0
for i in range(N-3):
#right boundary, diagonally above:
A[i*(N-1)+N-2][i*(N-1)+N-2+N] = 0
A = sp.csr_matrix(A)
return A
```
%% Cell type:code id: tags:
``` python
n = 2
N = pow(2,n)
h = 1/N
```
%% Cell type:code id: tags:
``` python
A = sparse_A(h).todense()
A
```
%% Output
matrix([[1024., -512., 0., -512., 256., 0., 0., 0., 0.],
[-512., 1024., -512., 256., -512., 256., 0., 0., 0.],
[ 0., -512., 1024., 0., 256., -512., 0., 0., 0.],
[-512., 0., 0., 1024., -512., 0., -512., 256., 0.],
[ 256., -512., 256., -512., 1024., -512., 256., -512., 256.],
[ 0., 256., -512., 0., -512., 1024., 0., 256., -512.],
[ 0., 0., 256., -512., 0., 0., 1024., -512., 0.],
[ 0., 0., 0., 256., -512., 256., -512., 1024., -512.],
[ 0., 0., 0., 0., 256., -512., 0., -512., 1024.]])
%% 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:code id: tags:
``` python
# for m inner nodes and l boundary nodes, B is a (m x l)-matrix
def sparse_B(h):
N = int(1/h)
m = pow(N-1,2)
l = 4*N
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, successor:
B[(i+1)*(N-1)-1][N + 2 * (i + 1)] = -2*pow(h,-4)
#upper boundary:
B[-(N-1)+i][-N+i] = -2*pow(h,-4)
#right above
B[-(N-1)+i][-N+i+1] = pow(h,-4)
#left above
B[-(N-2)+i][-N+i] = pow(h,-4)
for i in range(N-2):
#right boundary, diagonally above:
B[(i+1)*(N-1)-1][N + 2 * (i + 1)+2] = pow(h,-4)
#right boundary, diagonally below:
B[(i+2)*(N-1)-1][N + 2 * (i + 1)] = pow(h,-4)
B = sp.csr_matrix(B)
return B
```
%% Cell type:code id: tags:
``` python
# example for B
n = 3
h = pow(2,-n)
B = sparse_B(h)
B.todense()[-10:-1]
```
%% Output
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., 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., 4096., 0., -8192., 0., 4096., 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.,
-8192., 4096., 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.,
4096., -8192., 4096., 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., 4096., -8192., 4096., 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., 4096., -8192., 4096., 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., 4096., -8192., 4096., 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., 4096., -8192., 4096., 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):
N = int(1/h)
l = 4*N
g = np.zeros(l)
# 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
```
%% Cell type:code id: tags:
``` python
# example for g
n = 2
h = pow(2,-n)
g = vector_g(h)
g
```
%% 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:
``` python
n = 3
h = pow(2,-n)
N = pow(2,n)
```
%% Cell type:code id: tags:
``` python
A = sparse_A(h)
f = vector_f(h)
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:
``` python
RHS = f+B.dot(g)
appr_u = sp.linalg.spsolve(A,RHS)
```
%% Cell type:markdown id: tags:
#### **Error and Consistency**
%% 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
grid = np.zeros(7)
error_eukl = np.zeros(7)
error_inf = np.zeros(7)
for i in range(5):
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:
``` python
plt.loglog(grid,error_inf)
plt.loglog(grid, grid**2, label='2nd oder')
#plt.loglog(grid,grid, label='1st oder')
plt.ylabel('error in maximum-norm', fontsize=16)
plt.xlabel('step size h', fontsize=16)
```
%% Output
Text(0.5, 0, 'step size h')
%% Cell type:code id: tags:
``` python
plt.loglog(grid,error_eukl)
plt.loglog(grid, grid**2, label='2nd oder')
plt.loglog(grid,grid, label='1st oder')
plt.ylabel('error in Euklidean norm', fontsize=16)
plt.xlabel('step size h', fontsize=16)
```
%% Output
Text(0.5, 0, 'step size h')
%% Cell type:code id: tags:
``` python
```