问题:$$-\frac{d^2u}{dx^2}=f(x)$$
边值条件为$$u(a)=0,u(b)=1$$,
精确解为$$u=x^{\frac{3}{4}}$$
其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)=1
其双线性形式为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**(3/4)
def f_function(x):
"""
定义方程右端函数
"""
return 3/16*x**(-5/4)
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 hat_function(x,t,j,h,n):
#h=(b-a)/(n)
#t=[a+i*h for i in range(0,n+1)]
if j==0:
if x>=t[0] and x<=t[1]:
return 1-abs(x-t[0])/h[j]
else:
return 0
if j==n:
if t[n-1]<=x and x<=t[n]:
return 1-abs(x-t[n])/h[j]
else:
return 0
if j>0 and j<n:
if t[j-1]<=x and x<=t[j+1]:
return 1-abs(x-t[j])/h[j]
else:
return 0
"""
def compute_right_of_equation(point,weight,t,h,n,a,b):
"""
计算方程的右端项
参数:
h 步长列表
a,b 对应区间[a,b]
"""
c=(b-a)/2
s=(b-a)/2*point+(a+b)/2 #把区间[a,b]变化到[-1,1]
f=[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)]
f[-1]=f[-1]+1/h[-1] #添加右端初值条件
return f
def compute_stiffness_matrix(point,weight,t,h,n,a,b):
"""
计算系数矩阵,即总的刚度矩阵
"""
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 if x!=1]
error=function_u(np.array(t[:-1]))-uh
return uh,LA.norm(error,np.inf)
if __name__ == '__main__':
a=0 #区间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)+[1],'+b',label='approximate solution')
plt.plot(t,function_u(np.array(t)),'r',label='exact solution')
plt.title('Fig2.(where a=0,b=1,n=20,$u(x)=x^{3/4}$,$u(1)=0,u(1)=1$)')
plt.xlabel('x')
plt.ylabel('u')
plt.legend(loc='upper center')
print('all cost time:',time.time()-start_time)
运行结果
the error is: 4.12305198818e-05
all cost time: 0.23321843147277832