有限元求解两点边值问题之一

问题一：

$$\frac{d^2u}{dx^2}=f(x)$$

边界条件：$$u(a)=u(b)=0$$,

精确解为：$$u=x^4-1$$

运用Galerkin方法建立有限元方程(其过程可参考《微分方程数值解法》李荣华第四版，P226)

下面直接上代码：

# -*- coding: utf-8 -*-
"""
Created on Wed Jun  6 21:04:25 2018

@author: shaowu

Galerkin有限元法求解两边值问题，方程形式：
p*u''=f(u,x),u(a)=0,u(b)=0或u'(b)=0
其双线性形式为a(u,v)=p*u'*v'的一重积分
"""
import numpy as np
import pandas as pd
import scipy as sp
import numpy.linalg as LA
import matplotlib.pyplot as plt
from scipy.special.orthogonal import p_roots
import time
start_time=time.time()
def gauss_xw(m=100):
    """
    默认用100个点求高斯——勒让德节点point和权weight，并返回
    """
    point,weight=p_roots(m+1)
    return point,weight

def function_u(x):
    """精确解"""
    return x**4-1

def f_function(x):
    """
    定义方程右端函数
    """
    return 12*x**2

def hat_function(x,t,i,h,n):
    
    #构造试探函数空间的一组基函数：线性插值公式构造如下
    #参数：h 步长列表
    #t 节点列表
    #i 第i个节点
    #n 节点个数
    #x 变量
    if i==n:
        if t[n]>=x and x>=t[n-1]:
            return 1+(x-t[n])/h[n-1]
        else:
            return 0
    else:
        if t[i]>=x and x>=t[i-1]:
            return 1+(x-t[i])/h[i-1]
        elif t[i+1]>=x and x>=t[i]:
            return 1-(x-t[i])/h[i]
        else:
            return 0

def compute_right_of_equation(point,weight,t,h,n,a,b):
    """
    计算方程的右端项
    参数：
    h 步长列表
    a,b 对应区间[a,b]
    """
    #return [sp.integrate.quad(lambda x: f_function(x)*hat_function(x,t,j,h,n),t[j-1],t[j])[0]+\
    #        sp.integrate.quad(lambda x: f_function(x)*hat_function(x,t,j,h,n),t[j],t[j+1])[0] for \
    #        j in range(1,n)]+[sp.integrate.quad(lambda x: f_function(x)*hat_function(x,t,n,h,n),t[n-1],t[n])[0]]

    c=(b-a)/2
    s=(b-a)/2*point+(a+b)/2 #把区间[a,b]变化到[-1,1]
    return [sum([c*weight[i]*f_function(s[i])*hat_function(s[i],t,j,h,n) for i in range(len(s))])
            for j in range(1,len(t)-1)]
def compute_stiffness_matrix(point,weight,t,h,n,a,b):
    """
    计算系数矩阵，即总的刚度矩阵
    """
    #c=(b-a)/2
    #s=(b-a)/2*point+(a+b)/2 #把区间[a,b]变化到[-1,1]
    """
    A=[]
    for j in range(len(t)):
        A.append([sum([c*weight[k]*hat_function(s[k],t,i,h,n)*hat_function(s[k],t,j,h,n) for k in range(len(s))])
            for i in range(len(t))])
    """
    A=np.zeros((n-1,n-1))
    for j in range(1,len(t)-1):
        if j==1:
            A[j-1,j-1]=+sp.integrate.quad(lambda x: h[j]**(-2),t[j-1],t[j+1])[0]
                   #sp.integrate.quad(lambda x: h[j]**(-2),t[j],t[j+1])[0]
            A[j-1,j]=sp.integrate.quad(lambda x: -h[j]**(-2),t[j],t[j+1])[0]
            #A[j,j-1]=sp.integrate.quad(lambda x: h[j]**(-2),t[j-1],t[j])[0]
        
        elif j==n-1:
            A[j-1,j-2]=sp.integrate.quad(lambda x: -h[j-1]**(-2),t[j-1],t[j])[0]
            A[j-1,j-1]=+sp.integrate.quad(lambda x: h[j-1]**(-2),t[j-1],t[j+1])[0]
        else:
            A[j-1,j-2]=sp.integrate.quad(lambda x: -h[j-1]**(-2),t[j-1],t[j])[0]
            A[j-1,j-1]=+sp.integrate.quad(lambda x: h[j-1]**(-2),t[j-1],t[j])[0]+\
                   sp.integrate.quad(lambda x: h[j]**(-2),t[j],t[j+1])[0]
            A[j-1,j]=sp.integrate.quad(lambda x: -h[j]**(-2),t[j],t[j+1])[0]
    return np.array(A)

def solve_ui(A,b):
    """
    计算ui
    """
    return np.linalg.solve(A,-np.array(b))
def solve_uh(ui,t,h,n,a,b):
    """
    计算uh和误差
    """
    uh=[sum([ui[j]*hat_function(x,t,j+1,h,n) for j in range(n-1)]) for x in t]
    error=function_u(np.array(t))-uh
    #print(LA.norm(error,np.inf))
    return uh,LA.norm(error,np.inf)
if __name__ == '__main__':
    a=-1 #区间a,b
    b=1
    n=20
    h=(b-a)/(n) #等距步长
    t=[a+i*h for i in range(0,n+1)] #节点
    h=[(b-a)/(n) for i in range(n)] #子区间步长集合
    
    point,weight=gauss_xw(m=100) #高斯勒让德节点point和权weight
    f=compute_right_of_equation(point,weight,t,h,n,a,b) #计算方程右端项
    A=compute_stiffness_matrix(point,weight,t,h,n,a,b) #计算系数矩阵，即总刚度矩阵
    ui=solve_ui(A,f) #计算ui
    uh,error=solve_uh(ui,t,h,n,a,b) #计算uh和误差
    print('the error is:',error)
    plt.plot(t,[0]+list(ui)+[0],'+b',label='approximate solution')
    plt.plot(t,function_u(np.array(t)),'r',label='exact solution')
    plt.title('Fig1.(where a=-1,b=1,n=20,$u(x)=x^4-1$,$u(-1)=u(1)=0$)')
    plt.xlabel('x')
    plt.ylabel('u')
    plt.legend(loc='upper center')
    print('all cost time:',time.time()-start_time)

运行结果：

the error is: 0.000169677784248

all cost time: 0.3332369327545166

如有疑问，欢迎提问。。。