Python finite difference method for differential equations
Question: I must solve the Euler Bernoulli differential beam equation which is:

    
    
    u’’’’(x) = f(x) ; 
    

(x is the coordinate of the beam axis points)

and boundary conditions:

    
    
    u(0)=0,    u’(0)=0,    u’’(1)=0,    u’’’(1)=a 
    

I have studied the theory of numerically finite differences which expresses
the series of derivations as:

    
    
    U’k = (1/2*h)(Uk+1 - Uk-1) 
    
    U’’k = (1/h2)(Uk+1 - 2 Uk + Uk-1) 
    
    U’’’k = (1/2h3)(Uk+2 - 2 Uk+1 + 2 Uk-1 + Uk-2)
    
    U’’’’k = (1/h4)(Uk+2 - 4 Uk+1 + 6 Uk - 4 Uk-1 + Uk-2)
    

(`k+1`, `k+2`, etc. etc are subscripts)

and I found a script which expresses it as follows:

    
    
    import numpy as np
    from scipy.linalg import solveh_banded
    import matplotlib.pyplot as plt
    
    def beam4(n,ffun,a):
        x = np.linspace(0,1,n+1)
    
        h = 1.0/n
    
        stencil = np.array((1,-4,6))
    
        B = np.outer(stencil,np.ones(n))
    
        B[1,-1] = -2;    B[2,0] = 7;    B[2,-1] = 1;    B[2,-2] = 5
    
        f = ffun(x)
    
        f *= h** 4;     f[-1] *= 0.5;     f[-1] -= a*h**3
    
        u = np.zeros(n+1) 
    
        u[1:] = solveh_banded(B,f[1:],lower=False)
    
        return x,u
    

But I can't understand why the coefficient matrix is built this way:

    
    
    stencil = np.array((1,-4,6))
    
    B = np.outer(stencil,np.ones(n))
    
    B[1,-1] = -2;    B[2,0] = 7;    B[2,-1] = 1;    B[2,-2] = 5    
    
    f = ffun(x)
    
    f *= h**4;    f[-1] *= 0.5;    f[-1] -= a*h**3  "
    

Thanks in advance!!!!


Answer: Hope this is helpful. (due to this is my first time to post an answer)

The coefficients of this Hermitian positive-definite banded matrix are due to
applied of ghost node method. It is one of most efficient and popular method
for treating the boundary conditions of FDM without lossing of accuracy (here
these coefficients will give a second order converge rate in general). If you
have trouble to visual the matrix please check the 'K' matrix in my code
below:

    
    
    from numpy import linspace,zeros,array,eye,dot
    from numpy.linalg import solve
    from pylab import plot,xlabel,ylabel,legend,show
    a       = 0.2          ;b     = 0.0
    LX,dx   = 1.0,0.05     ;nx    = int(LX/dx)+1
    X       = linspace(0.0,LX,nx)
    Fs      = X**2.0
    """"""""""""""""""""""""""""""""""""""""""""""""""""""
    def calcB(l,b):
        if   b==0:    return [0.0,0.0]
        elif b==1:    return 1.0/l/l/l/l*array([-4.0,7.0,-4.0,1.0])
        elif b==nx-2: return 1.0/l/l/l/l*array([1.0,-4.0,5.0,-2.0])
        elif b==nx-1: return 1.0/l/l/l/l*array([2.0,-4.0,2.0])
        else:         return 1.0/l/l/l/l*array([1.0,-4.0,6.0,-4.0,1.0])
    U       = zeros(nx)     ;V       = zeros(nx)
    M       = eye(nx)       ;K       = zeros((nx,nx))
    F       = zeros(nx)     ;F[nx-2] = -b/dx/dx
    F[nx-1] = 2.0*b/dx/dx-2.0*a/dx
    for i in range(nx):
        if i == 0:      I   = [i,i+1]
        elif i == 1:    I   = [i-1,i,i+1,i+2]
        elif i == nx-2: I   = [i-2,i-1,i,i+1]
        elif i == nx-1: I   = [i-2,i-1,i]
        else:           I   = [i-2,i-1,i,i+1,i+2]
        for k,j in enumerate(I):
            K[i,j]  += calcB(dx,i)[k]
    """"""""""""""""""""""""""""""""""""""""""""""""""""""
    pn      = [0]        ;eq2     = pn
    eq1     = [i for i in range(nx) if i not in pn]
    
    MM1_    = K[eq1,:]; MM11,MM12 = MM1_[:,eq1],MM1_[:,eq2]
    RR      = F+Fs
    
    U[eq2]  = [0.0]
    U[eq1]  = solve(MM11,(RR[eq1]-dot(MM12,U[eq2]))) 
    ######################Plotting#########################
    Us = lambda x: x**6.0/360.0+x**3.0/6.0*(a-1.0/3.0)+x**2.0/2.0*(1.0/4.0-a)
    plot(X,U,'bo',label='FDM') ;plot(X,Us(X),'g-',label='solution')
    xlabel('X'); ylabel('U'); legend(loc='best')
    show()
    

cheers

