1-D_model

# -*- coding: utf-8 -*-
"""
fea method for a one-dimensional model problem:

require:
- numpy:    1.8.1
- scipy:    0.14.0

Yung Wong <yungwong<yungwong.seu@gmail.com>
May 26 09:26:38 2014
"""

import numpy as np
from scipy import integrate

h = 0.0001
M  = (1-0)/h-1
xi = np.linspace(0, 1, num=M+2)
f = lambda x: 4*np.pi**2*np.sin(2*np.pi*x)
A = (np.diag(2*np.ones(M))+np.diag(-1*np.ones(M-1),k=1)+\
    np.diag(-1*np.ones(M-1),k=-1))/h
fi = np.zeros(M)
for i in np.arange(1, M+1):
    def phi(x):
        if x < xi[i]:
            return 1/h*(x-xi[i-1])
        else:
            return -1/h*(xi[i+1]-x)
    fi[i-1], err = integrate.quad(lambda x: phi(x)*f(x), xi[i-1], xi[i+1])
p = np.linalg.solve(A, fi)
p = np.insert(p, 0, 0)
p = np.append(p, 0)