In [11]:
import scipy as sp
import numpy as np
import scipy.sparse as sparse
import scipy.sparse.linalg as sla
import matplotlib.pyplot as plt
%matplotlib inline
In [12]:
def hnorm(r):
    """define ||r||_h = h ||r||_2"""
    n = len(r)
    h = 1.0 / (n+1)
    hrnorm = h * np.linalg.norm(r)
    return hrnorm

def poissonop(n):
    A = (n+1)**2 * sparse.diags([-1, 2, -1], [-1, 0, 1], shape=(n,n), format='csr')
    return A

def residual(u, f):
    r = np.zeros(len(u))
    r[1:-1] = f[1:-1] - (n+1)**2 * (2 * u[1:-1] - u[2:] - u[:-2])
    r[0]    = f[0]    - (n+1)**2 * (2 * u[0]    - u[1])
    r[-1]   = f[-1]   - (n+1)**2 * (2 * u[-1]   - u[-2])
    return r

def relax(u, f, nu):
    n = len(u)
    Dinv = 1.0 / (2.0 * (n+1)**2)
    omega = 2.0 / 3.0
    unew = u.copy()
    
    for steps in range(nu):
        unew = unew + omega * Dinv * residual(unew, f)
    
    return unew

def interpolate(uc):
    """interpolate u of size 2**(k-1)-1 to 2**(k)-1"""
    uf = np.zeros((2*len(uc) + 1,))
    uf[:-1:2] = 0.5 * uc
    uf[1::2] = uc
    uf[2::2] += 0.5 * uc
    return uf

def restrict(uf):
    """restrict u of size 2**(k)-1 to 2**(k-1)-1"""
    uc = 0.25 * uf[:-1:2] + 0.5 * uf[1::2] + 0.25 * uf[2::2]
    return uc
In [19]:
def vcycle(u, f, nu):
    n = len(u)
    
    if n == 1:
        u = 0.125 * f
    elif n > 1:
        u = relax(u, f, nu)
        r = residual(u, f)
        rc = restrict(r)
        uc = np.zeros(rc.shape)
        
        #uc = vcycle(uc, rc, nu)
        Ac = poissonop(len(uc))
        uc = sparse.linalg.spsolve(Ac, rc)
        
        e = interpolate(uc)
        u = u + e        
        u = relax(u, f, nu)
      
    return u
In [22]:
k = 6
n = 2**k - 1
xx = np.linspace(0,1,n+2)[1:-1]

f = 2 - 12 * xx**2
ustar = xx**4 - xx**2

#f = np.pi**2 * np.sin(np.pi * xx)
#ustar = np.sin(np.pi * xx)
u = np.random.rand(len(f))

A = poissonop(n)
udiscretestar = sla.spsolve(A, f)

print(hnorm(residual(u, f)))

for i in range(25):
    u = vcycle(u, f, 1)
    res = hnorm(residual(u, f))
    err = hnorm(ustar - u)
    aerr = hnorm(udiscretestar - u)
    print("res = %10.4e, total err = %10.4e  alg err = %10.4e" % (res, err, aerr))

364.51259364
res = 2.3458e+01, total err = 5.0219e-03  alg err = 5.0202e-03
res = 2.6065e+00, total err = 5.5957e-04  alg err = 5.5780e-04
res = 2.8961e-01, total err = 6.3942e-05  alg err = 6.1977e-05
res = 3.2178e-02, total err = 1.0124e-05  alg err = 6.8864e-06
res = 3.5754e-03, total err = 5.8566e-06  alg err = 7.6515e-07
res = 3.9726e-04, total err = 5.5989e-06  alg err = 8.5017e-08
res = 4.4141e-05, total err = 5.5747e-06  alg err = 9.4463e-09
res = 4.9045e-06, total err = 5.5721e-06  alg err = 1.0496e-09
res = 5.4494e-07, total err = 5.5718e-06  alg err = 1.1662e-10
res = 6.0549e-08, total err = 5.5717e-06  alg err = 1.2958e-11
res = 6.7277e-09, total err = 5.5717e-06  alg err = 1.4397e-12
res = 7.4753e-10, total err = 5.5717e-06  alg err = 1.5993e-13
res = 8.3059e-11, total err = 5.5717e-06  alg err = 1.7775e-14
res = 9.2311e-12, total err = 5.5717e-06  alg err = 1.9492e-15
res = 1.0274e-12, total err = 5.5717e-06  alg err = 2.2424e-16
res = 1.1523e-13, total err = 5.5717e-06  alg err = 1.3085e-16
res = 1.4200e-14, total err = 5.5717e-06  alg err = 8.6159e-17
res = 1.0089e-14, total err = 5.5717e-06  alg err = 9.0040e-17
res = 3.8241e-15, total err = 5.5717e-06  alg err = 1.1134e-16
res = 3.4694e-18, total err = 5.5717e-06  alg err = 1.1319e-16
res = 0.0000e+00, total err = 5.5717e-06  alg err = 1.1319e-16
res = 0.0000e+00, total err = 5.5717e-06  alg err = 1.1319e-16
res = 0.0000e+00, total err = 5.5717e-06  alg err = 1.1319e-16
res = 0.0000e+00, total err = 5.5717e-06  alg err = 1.1319e-16
res = 0.0000e+00, total err = 5.5717e-06  alg err = 1.1319e-16
In [ ]: