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

error is still big

parent 13f39215
Branches
No related tags found
No related merge requests found
%% Cell type:markdown id: tags:
%% Cell type:code id: tags:
``` python
import numpy as np
```
%% Cell type:markdown id: tags:
First, the components of the following equation will be assembled:
$$A \underline{u} = \underline{f} + B\underline{g}$$
**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
def f(x,y):
#right-hand-side of the poisson equation
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))
```
%% Cell type:code id: tags:
%% Cell type:markdown id: tags:
``` python
n = 3
h = pow(2,-n)
N = pow(2,n)
```
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)))
return A
```
%% Cell type:code id: tags:
``` python
matrix_A(h).shape
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
```
%% Output
(49, 49)
%% Cell type:code id: tags:
``` python
# 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:code id: tags:
``` python
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:
``` python
n = 3
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)
```
%% 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
u= 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:
``` python
u[-1]
```
%% Output
-13.228477586622512
%% Cell type:code id: tags:
``` python
appr_u[-1]
```
%% Output
0.025380345287435015
%% Cell type:code id: tags:
``` python
f.shape
```
%% Output
(16129,)
%% Cell type:code id: tags:
``` python
A.shape
```
%% Output
(16129, 16129)
%% Cell type:code id: tags:
``` python
u.shape
```
%% Output
(16129,)
%% Cell type:code id: tags:
``` python
ex=np.matmul(A,u)
```
%% Cell type:code id: tags:
``` python
np.eye?
ex[-1]
```
%% Output
434807.0641872273
%% Cell type:code id: tags:
``` python
```
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment