In [1]:
import numpy as np
import matplotlib.pyplot as plt
import scipy as sp
import scipy.sparse as sparse
%matplotlib inline
In [68]:
def uex(x):
    return np.sin(np.pi * x / 2)

def f(x):
    return (np.pi / 2)**2 * np.sin(np.pi * x / 2)
In [79]:
n = 20
x = np.linspace(0,1,n+1)
dx = x[1] - x[0]
I = []
J = []
AA = []
If = []
Jf = []
FF = []
for i in range(0,n):
    row = [i, i, i+1, i+1]
    col = [i, i+1, i, i+1]
    data = [1, -1, -1, 1]
    I += row
    J += col
    AA += data
    
    row = [i, i+1]
    col = [0, 0]
    data = [f(x[i]), f(x[i+1])]
    If += row
    Jf += col
    FF += data
    
A = (1/dx) * sparse.coo_matrix((AA, (I, J)))
F = (dx/2) * sparse.coo_matrix((FF, (If, Jf))).todense()
F = np.array(F).ravel()
In [80]:
plt.plot(x,F, x, dx* f(x))
F
Out[80]:
array([ 0.        ,  0.0096795 ,  0.01929933,  0.02880017,  0.03812344,
        0.04721168,  0.05600883,  0.06446068,  0.0725151 ,  0.08012244,
        0.0872358 ,  0.09381133,  0.09980847,  0.10519026,  0.10992352,
        0.11397907,  0.11733189,  0.11996133,  0.12185116,  0.12298975,
        0.06168503])
In [81]:
A = A.tocsr()
A[0,0] = 1
A[0,1] = 0
A[1,0] = 0
F[0] = 0
A.todense()
Out[81]:
matrix([[  1.,   0.,   0.,   0.,   0.,   0.,   0.,   0.,   0.,   0.,   0.,
           0.,   0.,   0.,   0.,   0.,   0.,   0.,   0.,   0.,   0.],
        [  0.,  40., -20.,   0.,   0.,   0.,   0.,   0.,   0.,   0.,   0.,
           0.,   0.,   0.,   0.,   0.,   0.,   0.,   0.,   0.,   0.],
        [  0., -20.,  40., -20.,   0.,   0.,   0.,   0.,   0.,   0.,   0.,
           0.,   0.,   0.,   0.,   0.,   0.,   0.,   0.,   0.,   0.],
        [  0.,   0., -20.,  40., -20.,   0.,   0.,   0.,   0.,   0.,   0.,
           0.,   0.,   0.,   0.,   0.,   0.,   0.,   0.,   0.,   0.],
        [  0.,   0.,   0., -20.,  40., -20.,   0.,   0.,   0.,   0.,   0.,
           0.,   0.,   0.,   0.,   0.,   0.,   0.,   0.,   0.,   0.],
        [  0.,   0.,   0.,   0., -20.,  40., -20.,   0.,   0.,   0.,   0.,
           0.,   0.,   0.,   0.,   0.,   0.,   0.,   0.,   0.,   0.],
        [  0.,   0.,   0.,   0.,   0., -20.,  40., -20.,   0.,   0.,   0.,
           0.,   0.,   0.,   0.,   0.,   0.,   0.,   0.,   0.,   0.],
        [  0.,   0.,   0.,   0.,   0.,   0., -20.,  40., -20.,   0.,   0.,
           0.,   0.,   0.,   0.,   0.,   0.,   0.,   0.,   0.,   0.],
        [  0.,   0.,   0.,   0.,   0.,   0.,   0., -20.,  40., -20.,   0.,
           0.,   0.,   0.,   0.,   0.,   0.,   0.,   0.,   0.,   0.],
        [  0.,   0.,   0.,   0.,   0.,   0.,   0.,   0., -20.,  40., -20.,
           0.,   0.,   0.,   0.,   0.,   0.,   0.,   0.,   0.,   0.],
        [  0.,   0.,   0.,   0.,   0.,   0.,   0.,   0.,   0., -20.,  40.,
         -20.,   0.,   0.,   0.,   0.,   0.,   0.,   0.,   0.,   0.],
        [  0.,   0.,   0.,   0.,   0.,   0.,   0.,   0.,   0.,   0., -20.,
          40., -20.,   0.,   0.,   0.,   0.,   0.,   0.,   0.,   0.],
        [  0.,   0.,   0.,   0.,   0.,   0.,   0.,   0.,   0.,   0.,   0.,
         -20.,  40., -20.,   0.,   0.,   0.,   0.,   0.,   0.,   0.],
        [  0.,   0.,   0.,   0.,   0.,   0.,   0.,   0.,   0.,   0.,   0.,
           0., -20.,  40., -20.,   0.,   0.,   0.,   0.,   0.,   0.],
        [  0.,   0.,   0.,   0.,   0.,   0.,   0.,   0.,   0.,   0.,   0.,
           0.,   0., -20.,  40., -20.,   0.,   0.,   0.,   0.,   0.],
        [  0.,   0.,   0.,   0.,   0.,   0.,   0.,   0.,   0.,   0.,   0.,
           0.,   0.,   0., -20.,  40., -20.,   0.,   0.,   0.,   0.],
        [  0.,   0.,   0.,   0.,   0.,   0.,   0.,   0.,   0.,   0.,   0.,
           0.,   0.,   0.,   0., -20.,  40., -20.,   0.,   0.,   0.],
        [  0.,   0.,   0.,   0.,   0.,   0.,   0.,   0.,   0.,   0.,   0.,
           0.,   0.,   0.,   0.,   0., -20.,  40., -20.,   0.,   0.],
        [  0.,   0.,   0.,   0.,   0.,   0.,   0.,   0.,   0.,   0.,   0.,
           0.,   0.,   0.,   0.,   0.,   0., -20.,  40., -20.,   0.],
        [  0.,   0.,   0.,   0.,   0.,   0.,   0.,   0.,   0.,   0.,   0.,
           0.,   0.,   0.,   0.,   0.,   0.,   0., -20.,  40., -20.],
        [  0.,   0.,   0.,   0.,   0.,   0.,   0.,   0.,   0.,   0.,   0.,
           0.,   0.,   0.,   0.,   0.,   0.,   0.,   0., -20.,  20.]])
In [82]:
import scipy.sparse.linalg as sla
In [83]:
u = sla.spsolve(A, F)
In [84]:
u
Out[84]:
array([ 0.        ,  0.07849944,  0.1565149 ,  0.2335654 ,  0.30917589,
        0.38288021,  0.45422394,  0.52276723,  0.58808749,  0.64978199,
        0.70747038,  0.76079697,  0.80943299,  0.85307859,  0.89146468,
        0.92435459,  0.95154555,  0.97286991,  0.98819621,  0.99742995,
        1.0005142 ])
In [85]:
plt.plot(x, np.abs(u - uex(x)))
Out[85]:
[<matplotlib.lines.Line2D at 0x10b1bbef0>]
In [ ]: