有限元编程：菜鸟篇
《有限元编程：菜鸟篇》
相信很多做过有限差分之后又想做做有限元的初学者会有和我一样的困惑，能看懂有限元算法的理论分析，但是真正应用到实际编程当中之前心里发怵，废话不多说，求人不如求己，看懂这篇文章将会让你迅速掌握有限元最基础的编程思想。
二、以经典扩散方程为例(反常扩散方程可类比此例)
考虑如下扩散方程初边值问题
\begin{equation*}\label{eq1.1}
\left \{
\begin{array}{ll}
\frac{\partial u(x,t)}{\partial t}=\frac{\partial^2u(x,t)}{\partial x^2}+f(x,t), \quad & a< x< b,0< t\leq T, \\
u(a,t)=0,\quad u(b,t)=0,\quad & 0\leq t\leq T, \\
u(x,0)=\varphi (x),\quad & a\leq x \leq b, \\
\end{array}
\right.\tag{1}
\end{equation*}
首先我们需要对时间与空间区间进行剖分，其中$M,N$分别记为空间与时间方向的剖分，$h=\frac{b-a}{M}$, $ \tau = \frac{T}{N} $，$h,\tau$分别为空间与时间方向的步长，则$a=x_0<x_1<\cdots < x_M=b$, $0=t_0<t_1<\cdots< t_N=T$。
空间方向利用有限元逼近， 可得（1）式的变分形式（弱形式）
$$\big( \frac{\partial u_h(t)}{\partial t},v_h \big)+a\big(u_h,v_h\big)=\big(f,v_h\big),\qquad v_n \in X_h \tag{2}$$.
其中时间方向利用向后差分离散,可得求解问题（1）的全离散格式
$$ \big( u^n_h,v_h \big)+\tau a\big(u^n_h,v_h\big)=\big( u^{n-1}_h,v_h \big)+\tau\big(f(x,t_n),v_h\big) \tag{3}$$
其中, 记双线形式$a(u^n_h,v_h)=\int_\Omega (u^n_h)'(v_h)'dx$为刚度矩阵（stiffness matrix），$(u^n_h,v_h)=\int_\Omega u^n_h v_h dx$为质量矩阵（mass matrix），$(f(x,t_n),v_h)=\int_\Omega f(x,t_n) v_hdx$ 为载荷向量(load vector)。$\Omega=[a,b]=\underset{i=1,2,3,\cdots M-1}{\cup} \Omega_i$，$X_h=\{ \phi_i \}$为有限元空间，这里我们选用的是线性插值
$$\phi_i=\left \{
\begin{array}{rl}
\frac{x-x_{i-1}}{h},&\qquad \mbox{当}\quad x\in \Omega_i,\\
\frac{x_{i+1}~-~x}{h},&\qquad \mbox{当}\quad x\in \Omega_{i+1},\\
0, &\qquad \mbox{其他}.\\
\end{array}
\right.\tag{4}$$
因此，我们需要获得的问题（1）的有限元解$u^n_h$可表示为
$$u^n_h=\sum\limits_{j=1}^{M-1} u_j^n \phi_j ,\tag{5}$$
其中$u^n_i$ 为问题（1）在点$(x_i,t_n)$处的数值解,将（5）式代入（3）式，我们可以获得求解问题（1）的有限元格式的代数方程组。.
$$\Big( \tau A + B \Big)u^n=Bu^{n-1}+\tau \widehat{f}^n.\tag{6}$$
其中$u^n=[u_1^n,u_2^n,\cdots,u_{M-1}^n]^T$, $A=(a_{ij})_{M-1\times M-1}$为刚度矩阵，由（4）式我们可以发现，对任意$i$, $\phi_i$只在$\Omega_i=[x_{i-1},x_i]$与$\Omega_{i+1}=[x_i,x_{i+1}]$两个区间不为零，其余区间即$|i-j|\geq 2$时$a_{ij}$均为零，则通过区间转换计算，我们发现
$$\begin{array}{rl}a_{i,i}&=\int_\Omega (\phi_i)'(\phi_i)'dx=\sum\limits_{i=1}^{M} \int_{\Omega_i} (\phi_i)'(\phi_i)'dx\\&=\int_{x_{i-1}}^{x_i} (\phi_i)'(\phi_i)'dx+\int_{x_i}^{x_{i+1}} (\phi_i)'(\phi_i)'dx\\&\underset{\Longrightarrow}{\mbox{区间转换}}\int_0^1 \frac{1}{h}d\xi+\int_0^1 \frac{1}{h} d\xi\\&=\frac{2}{h},\qquad i=1,2,\cdots, M-1.\tag{7}\end{array}$$
$$\begin{array}{rl}a_{i-1,i}=a_{i,i-1}&=\int_\Omega (\phi_i)'(\phi_{i-1})'dx=\sum\limits_{i=1}^{M} \int_{\Omega_i} (\phi_i)'(\phi_{i-1})'dx\\&=\int_{x_{i-1}}^{x_{i-1}} (\phi_i)'(\phi_i)'dx\\&=\int_0^1 \frac{-1}{h}d\xi\\&=\frac{-1}{h},\qquad i=2,3,\cdots,M-1.\tag{8}\end{array}$$
可知刚度矩阵$A$为一三对角矩阵，同理可知质量矩阵$B=(b_{ij})_{M-1\times M-1}$也为三对角矩阵，其矩阵元素由如下计算获得
$$\begin{array}{rl}b_{i,i}&=\int_\Omega (\phi_i)^2dx=\sum\limits_{i=1}^{M} \int_{\Omega_i} (\phi_i)^2dx\\&=\int_{x_{i-1}}^{x_i} (\phi_i)^2dx+\int_{x_i}^{x_{i+1}} (\phi_i)^2dx\\&=h \int_0^1 \xi^2d\xi+h\int_0^1 (1-\xi)^2d\xi\\&=\frac{2h}{3},\qquad i=1,2,\cdots, M-1.\tag{9}\end{array}$$
$$\begin{array}{rl}b_{i-1,i}=b_{i,i-1}&=\int_\Omega (\phi_i)(\phi_{i-1})dx=\sum\limits_{i=1}^{M} \int_{\Omega_i} (\phi_i)(\phi_{i-1})dx\\&=\int_{x_{i-1}}^{x_{i-1}} (\phi_i)(\phi_i)dx\\&=h\int_0^1 \xi(1-\xi)d\xi \\&=\frac{h}{6},\qquad i=2,3,\cdots,M-1.\tag{10}\end{array}$$
到了这里，我们终于可以松一口气了，主要的东西我们已经准备好了，，，，等等，，，还有载荷向量$\widehat{f}^n=[(f(x,t_n),\phi_1),(f(x,t_n),\phi_2),\cdots,(f(x,t_n),\phi_{M-1})]^T$还没计算，接下来给出载荷向量的计算
$$\begin{array}\Big( f(x,t_n),\phi_i \Big)&=\int_\Omega f(x,t_n)\phi_i dx=\sum\limits^{M-1}_{i=1} \int_{\Omega_{i}}f(x,t_n)\phi_i dx\\&=\int_{x_{i-1}}^{x_i} f(x,t_n)\phi_idx+\int_{x_i}^{x_{i+1}} f(x,t_n)\phi_idx\\&=h\int_0^1 f(x_{i-1}+h\xi,t_n)xdx+h\int_0^1 f(x_i+h\xi,t_n)(1-\xi)d\xi.\tag{11}\end{array}$$
为计算以上的一重积分，我在程序中采用的是17点$Gauss$求积公式，绝对符合我们对算法精度的要求！不会$Gauss$积分公式的同学可以查阅文献（黄云清等，数值计算方法，科学出版社，p.116）.具体程序实现如下
function u = single_quad( a,b,f )
% This rountine is written for Stiff/mass matrix terms, in which five-point
% Gauss numerical integral is used.
% f(x) is standar for integrand in interval [a,b]
% g(x) is transformed from f(x), and [a,b] turn into [-1,1].
M = 17;
g = gauss_transform( f,a,b );
y(1)=0.1784841814958478558506775;y(2)=0.3512317634538763152971855;y(3)=0.5126905370864769678862466;
y(4)=0.6576711592166907658503022;y(5)=0.7815140038968014069252301;y(6)=0.8802391537269859021229557;
y(7)=0.9506755217687677612227170;y(8)=0.9905754753144173356754340;y(9)=0;y(10)=-y(8);y(11)=-y(7);
y(12)=-y(6);y(13)=-y(5);y(14)=-y(4);y(15)=-y(3);y(16)=-y(2);y(17)=-y(1);
A(1)=0.1765627053669926463252710;A(2)=0.1680041021564500445099707;A(3)= 0.1540457610768102880814316;
A(4)=0.1351363684685254732863200;A(5)=0.1118838471934039710947884;A(6)=0.0850361483171791808835354;
A(7)=0.0554595293739872011294402;A(8)=0.0241483028685479319601100;A(9)= 0.1794464703562065254582656;
A(10)=A(8);A(11)=A(7);A(12)=A(6);A(13)=A(5);A(14)=A(4);A(15)=A(3);A(16)=A(2);A(17)=A(1);
u = 0;
for i = 1:M
u = u + A(i) * g( y(i) );
end
function g = gauss_transform( f,a,b )
g = @(y) ( b - a ) / 2 * f( ( a+b ) / 2 + ( b-a ) / 2 * y );
end
end
到了这里，我们基本清楚了有限元格式基本算法的实现,以下我逐一介绍代码的实现过程。
在编程之前，我们自己需要构思一下，我们到底需要干些什么事情，第一步：输入算例的基本参数，第二步：区间剖分，第三步全离散格式（5）的封装，第四步：数据处理（收敛阶及误差），可视化。要想比较有条理的一步一步实现，我建议将2-4步单独封装成一个函数文件，第一步写成脚本文件，因为第一步需要我们手动输入模型参数，当你换另一个算例时只需要在第一步的文件中修改参数，接下来看看第一步:(取真解$u(x,t)=x^2(x-2)e^t$, $f(x,t)=(x^3-2x^2-6x+4)e^t$, $[a,b]=[0,2], T=1$, 初值$u(x,0)=x^2(x-2)$.)
脚本文件：run_main.m
clear;clc
% author: Dongdong Hu
% date: 16-Sep-2018
% version:8.3.0.532 (R2014a)
M = 100; N =100; Lx = 0; Rx = 2; Lt = 0; Rt = 1;
left_condation = @( t ) 0;
right_condation = @( t ) 0;
initial_condation = @( x ) x.^2 .* ( x-2 );
f = @( x,t ) exp( t ) .* ( x.^3 - 2 * x.^2 - 6 * x + 4 );
exact = @( x,t ) exp( t ) .* x.^2 .* ( x-2 );
show = show_solution( );
LFEM = model_data( Lx,Rx,Lt,Rt,left_condation,right_condation,...
initial_condation,exact);
[ x,t,u ] = Linear_FEM( M,N,LFEM,f );
ue = LFEM.exact( x,t );
show.image_3D( x,t,u );
show.image_2D( x,u,ue )
show.image_2D_error( x,abs( u-ue ) )
max_error = show.max_error( [ 10,20,40,80 ],[ 10,20,40,80 ].^2,LFEM,f );
max_rate = show.error_rate( max_error )
[ L2_error,space_step ] = show.L2_error( [ 10,20,40,80 ],[ 10,20,40,80 ].^2,LFEM,f );
L2_rate = show.error_rate( L2_error )
show.plot_rate( max_error,space_step );
第二步，我们需要进行网格剖分及初边值条件等已知参数的输入，这些都属于模型的固有属性
函数文件：model_data.m
function LFEM = model_data( Lx,Rx,Lt,Rt,left_condation,right_condation,...
initial_condation,EXACT)
%MODEL_DATA 此处显示有关此函数的摘要
% 此处显示详细说明
LFEM = struct('space_grid',@space_grid,...
'time_grid',@time_grid,'left_boundary',@left_boundary,...
'right_boundary',@right_boundary,'initial_value',...
@initial_value,'exact',@exact);
function [ x,h ] = space_grid( M )
h = ( Rx - Lx ) / M;
x = linspace( Lx,Rx,M+1 );
end
function [t,tau] = time_grid( N )
tau = ( Rt - Lt ) / N;
t = linspace( Lt,Rt,N+1 );
end
function u = left_boundary( t )
u = left_condation( t );
end
function u = right_boundary( t )
u = right_condation( t );
end
function u = initial_value( x )
u = initial_condation( x );
end
function u = exact( x,t )
u = zeros( length( x ),length( t ) );
for k = 1:length( t )
u( 1:end,k ) = EXACT( x,t( k ) );
end
end
end
第三步，我们需要将算法（5）封装
函数文件：Linear_FEM.m
function [ x,t,u ] = Linear_FEM( M,N,LFEM,f )
%L1_LINEAR_FEM Galerkin Finite elemeth method
% this programming apply linear interpolation in finite elment method
% with respect to x, the integer-order time partial derivative is
% approximated by backward difference method.
u = zeros( M+1,N+1 );
[ x,h ] = LFEM.space_grid( M );
[t,tau] = LFEM.time_grid( N );
u( 1,1:N+1 ) = LFEM.left_boundary( t );
u( M+1,1:N+1 ) = LFEM.right_boundary( t );
u( 1:M+1,1 ) = LFEM.initial_value( x );
stiffness_matrix = zeros( M-1 );
mass_matrix = zeros( M-1 );
for i = 1:M-1
stiffness_matrix( i,i ) = 2 / h;
mass_matrix( i,i ) = 2 / 3 * h;
if i>1
stiffness_matrix(i-1,i) = -1 / h;
mass_matrix(i-1,i) = h / 6;
stiffness_matrix(i,i-1) = -1 / h;
mass_matrix(i,i-1) = h / 6;
end
end
for n = 1:N
nn = n + 1;
load_vector = tau * h * single_quad( 0,1,@( xi ) f( x( 1:end-2 ) + h * xi,t( nn ) ) .* xi )'...
+ tau * h * single_quad( 0,1,@( xi ) f( x( 2:end-1 ) + h * xi,t( nn ) ) .* ( 1-xi ) )'...
+ mass_matrix * u( 2:end-1,nn-1 );
u( 2:end-1,nn ) = ( tau * stiffness_matrix + mass_matrix ) \ ( load_vector );
end
end
第四步，我们需要对我们所获得的数值解进行数据分析
函数文件：show_solution.m
function show = show_solution( )
%SHOW_SOLUTION 此处显示有关此函数的摘要
% 此处显示详细说明
show = struct( 'image_3D',@image_3D,'image_2D',@image_2D,'image_2D_error',@image_2D_error,...
'max_error',@max_error,'error_rate',@error_rate,'L2_error',@L2_error,'L2_rate',@L2_rate,...
'plot_rate',@plot_rate );
function image_3D( X,T,u )
figure
[ x,t ] = meshgrid( T,X );
mesh( x,t,u );
xlabel('x');ylabel('y');zlabel('u');
end
function image_2D( x,u,ue )
figure
plot( x,u(1:end,end),'b*' );hold on;
plot( x,ue( 1:end,end ),'r-' );
legend('numerical.sol','exact.sol');
xlabel('x');ylabel('u( x,t=1 )');
end
function image_2D_error( x,u )
figure
plot( x,u(1:end,end) );
xlabel('x');ylabel('|error( x,t=1 )|');
end
function Max_error = max_error( M,N,LFEM,f )
Max_error = zeros( size( M ) );
for i = 1:length( M )
[ x,t,u ] = Linear_FEM( M(i),N(i),LFEM,f );
ue = LFEM.exact( x,t );
Max_error( i ) = max( max( abs( u-ue ) ) );
end
end
function Error_rate = error_rate( max_error )
Error_rate = log2( max_error(1:end-1) ./ max_error( 2:end ) );
end
function [ L2_Error,space_step ] = L2_error( M,N,LFEM,f )
L2_Error = zeros( size( M ) );
space_step = L2_Error;
for i = 1:length( M )
[ x,t,u ] = Linear_FEM( M(i),N(i),LFEM,f );
ue = LFEM.exact( x,t );
[ ~,h ] = LFEM.space_grid( M(i) );
[~,tau] = LFEM.time_grid( N(i) );
space_step( i ) = h;
L2_Error(i) = sqrt( h * tau * sum( sum( ( u( 2:end-1, 2:end )-ue( 2:end-1, 2:end ) ).^2 ) ) );
end
end
function plot_rate( error1,error2,space_step )
figure
loglog( space_step,error1,'-kp' );hold on;
loglog( space_step,error2,'-r^' );hold on;
loglog( space_step,space_step.^2,'-g*' );
legend('||error||_{\infty}','|| error ||_{L_2}','Ch^2');
xlabel('log(1/h)'); ylabel('error');
end
end
到这里，整个编程到此结束，下面我来展示程序运行结果
请注意，虽然数值结果看似比较合理，但是，我在程序中计算$|| \cdot ||_{L_2}$范数意义下的误差是犯了一点错误的，我们知道在有限差分当中计算$|| \cdot ||_{L_2}$范数误差时用的是离散形式的定义，但是在有限元当中我们应该严格使用连续的$|| \cdot ||_{L_2}$范数定义：
$$\begin{array}{rl} || u(x,t_n)-u^n_h||_{L_2}&=\int_{\Omega} (u^n-u^n_h)^2dx\\&=\sum\limits^{M}_{i=1} \int_{\Omega_i} (u^n-u_{h,i}^n)^2dx \\&=\sum\limits^{M}_{i=1} \int_{x_{i-1}}^{x_i} (u^n-u_{h,i}^n)^2dx\\&=\sum\limits^{M}_{i=1} h\int_0^1 \Big(u^n( x_{i-1}+h\xi )-u_{i-1}^n(1-\xi)-\xi u_i^n \Big)^2d\xi. \tag{12}\end{array}$$
由（5）式我们可得
$$ u_{h,i}^n(x)=\frac{x_i-x}{h}u_{i-1}^n+\frac{x-x_{i-1}}{h}u_i^n \qquad x\in \Omega_i$$
通常我们取$n=N$, 也就是测量最后一时刻的观测阶.（综上分析，我们可以将第四部文件及相应处代码做一点修改，请参考代码：点击下载>>）
以上仅仅是是最简单的有限元编程（有不当之处，望留言指正），下一篇将不定期更新$caputo$型的时间分数阶扩散方程的有限元编程，一维空间分数阶扩散方程的有限元编程，二维空间分数阶扩散方程的有限元编程，二维时-空分数阶扩散方程的有限元编程，二维时-空分数阶扩散方程的有限元ADI算法编程
