Skip to content
Snippets Groups Projects
Commit 9c6c69bf authored by penrose's avatar penrose
Browse files

added boundary to matrix A, still an error at the corner of 1,1

parent b7757101
No related branches found
No related tags found
No related merge requests found
%% 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: %% 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: %% Cell type:code id: tags:
``` python ``` python
import numpy as np #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: %% 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}$$
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
def u(x,y): def matrix_A(h):
return pow(x,4)*pow(y,5)-17*np.sin(x*y) 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: %% Cell type:code id: tags:
``` python ``` python
def f(x,y): def vector_f(h):
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)) 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:
# just the initialisation of matrix B
def matrix_B(h):
N = int(1/h)
m = pow(N-1,2)
l = 4*N
B = np.zeros((m,l))
return B
%% Cell type:raw id: tags:
def vector_g(h):
N = int(1/h)
l = 4*N
g = np.zeros(l)
g[-1] = u(1,1)
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: %% Cell type:code id: tags:
``` python ``` python
n = 2 n = 6
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
def matrix_A(h): A=matrix_A(h)
N = int(1/h) f=vector_f(h)
m = pow(N-1,2) appr_u = np.linalg.solve(A,f)
A = np.zeros((m,m))
return A
``` ```
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
matrix_A(h) A[3968,3905]
``` ```
%% Output %% Output
array([[0., 0., 0., 0., 0., 0., 0., 0., 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., 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.]])
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
def matrix_B(h): def exact_solution(h):
N = int(1/h) N = int(1/h)
m = pow(N-1,2) l = pow(N-1,2)
l = 4*N v = np.zeros(l)
B = np.zeros((m,l)) for i in range(N-1):
return B 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: %% Cell type:code id: tags:
``` python ``` python
def vector_g(h): v= exact_solution(h)
N = int(1/h) ```
l = 4*N
g = np.zeros(l) %% Cell type:code id: tags:
g[-1] = u(1,1)
return g ``` python
error = max(abs(appr_u - v))
```
%% Cell type:code id: tags:
``` python
appr_u[-1]
```
%% Output
0.015864236841443172
%% Cell type:code id: tags:
``` python
error
```
%% Output
13.161396848144852
%% Cell type:code id: tags:
``` python
``` ```
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
``` ```
......
%% 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
from scipy.sparse.linalg import spsolve
``` ```
%% 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 is given 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 \; \partial \Omega$$
%% 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 -(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))
``` ```
%% 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}$$
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
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)))
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 return A
``` ```
%% 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:raw id: tags:
``` python
# just the initialisation of matrix B # just the initialisation of matrix B
def matrix_B(h): def matrix_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))
return B return B
```
%% Cell type:code id: tags: %% Cell type:raw id: tags:
``` python
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) g[-1] = u(1,1)
return g return g
```
%% Cell type:markdown id: tags: %% Cell type:markdown id: tags:
Since the exact solution is zero at the boundary, the product B*g is zero. Since the exact solution is zero at the boundary, the product B*g is zero.
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
n = 3 n = 6
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=matrix_A(h)
f=vector_f(h) f=vector_f(h)
appr_u = np.linalg.solve(A,f) appr_u = np.linalg.solve(A,f)
``` ```
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
A[3968,3905]
```
%% Output
4096.0
%% Cell type:code id: tags:
``` 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
u= exact_solution(h) v= exact_solution(h)
``` ```
%% Output
---------------------------------------------------------------------------
TypeError Traceback (most recent call last)
<ipython-input-110-a5ff9eca95bf> in <module>
----> 1 u= exact_solution(h)
<ipython-input-109-71aecd4cad10> in exact_solution(h)
5 for i in range(N-1):
6 for k in range(N-1):
----> 7 v[k+i*(N-1)]=u((k+1)/(N),(i+1)/(N))
8 return v
TypeError: 'numpy.ndarray' object is not callable
%% Cell type:code id: tags:
``` python
max(abs(appr_u - u))
```
%% Output
---------------------------------------------------------------------------
ValueError Traceback (most recent call last)
<ipython-input-106-a6b4164c5b68> in <module>
----> 1 max(abs(appr_u - u))
ValueError: operands could not be broadcast together with shapes (49,) (16129,)
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
u[-1] error = max(abs(appr_u - v))
``` ```
%% Output
-13.228477586622512
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
appr_u[-1] appr_u[-1]
``` ```
%% Output %% Output
0.025380345287435015 0.015864236841443172
%% Cell type:code id: tags:
``` python
f.shape
```
%% Output
(16129,)
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
A.shape error
``` ```
%% Output %% Output
(16129, 16129) 13.161396848144852
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
u.shape
``` ```
%% Output
(16129,)
%% Cell type:code id: tags:
``` python
ex=np.matmul(A,u)
```
%% Cell type:code id: tags:
``` python
ex[-1]
```
%% Output
434807.0641872273
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
``` ```
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment