《有限元編程：菜鳥篇》
相信很多做過有限差分之后又想做做有限元的初學者會有和我一樣的困惑，能看懂有限元算法的理論分析，但是真正應用到實際編程當中之前心里發怵，請教學過有限元程序的同學的時候，他們往往會，這個怎么怎么的簡單，這個你怎么能不會？這個不就是什么什么嗎bulabula...這時候你的心里一定和我一樣有一萬匹草泥馬在心里奔騰，廢話不多說，求人不如求己，這篇文章將會讓你迅速掌握有限元最基礎的編程思想。
二、以經典擴散方程為例(反常擴散方程可類比此例)
考慮如下擴散方程初邊值問題
\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算法編程