快速学会有限元编程
Posted xtu-hudongdong
tags:
篇首语：本文由小常识网(cha138.com)小编为大家整理，主要介绍了快速学会有限元编程相关的知识，希望对你有一定的参考价值。
《快速学会有限元编程》
一、前言
相信很多做过有限差分之后又想做做有限元的初学者会有和我一样的困惑，能看懂有限元算法的理论分析，但是真正应用到实际编程当中之前心里发怵，请教学过有限元程序的同学的时候，他们往往会，这个怎么怎么的简单，这个你怎么能不会？这个不就是什么什么吗bulabula...这时候你的心里一定和我一样有一万匹草泥马在心里奔腾，我特么的要是会还来问你，废话不多说，求人不如求己，这篇文章将会让你迅速掌握有限元最基础的编程思想。
二、以经典扩散方程为例
考虑如下扩散方程初边值问题
egin{equation*}label{eq1.1}
left {
egin{array}{rl}
frac{partial u(x,t)}{partial t}=frac{partial^2u(x,t)}{partial t}+f(x,t), quad & a< x< b,0< tleq T, \\
u(a,t)=0,quad u(b,t)=0,quad & 0leq tleq T, \\
u(x,0)=varphi (x),quad & aleq x leq b, \\
end{array}
ight. ag{1}
end{equation*}
首先我们需要对时间与空间区间进行剖分，其中$M,N$分别记为空间与时间方向的剖分，$h=frac{b-a}{M}$, $ au = frac{T}{N} $，$h, au$分别为空间与时间方向的步长，则$a=x_0<x_1<cdots < x_M=b$, $0=t_0<t_1<cdots< t_N=T$。
空间方向利用有限元逼近， 可得（1）式的变分形式（弱形式）
$$ig( frac{partial u_h(t)}{partial t},v_h ig)+aig(u_h,v_hig)=ig(f,v_hig),qquad v_n in X_h ag{2}$$.
其中时间方向利用向后差分离散,可得求解问题（1）的全离散格式
$$ ig( u^n_h,v_h ig)+ au aig(u^n_h,v_hig)=ig( u^{n-1}_h,v_h ig)+ auig(f^n,v_hig) ag{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^n,v_h)=int_Omega f^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 {
egin{array}{ll}
frac{x_{i+1}-x}{h},qquad xin Omega_i,\\
frac{x-x_{i-1}}{h},qquad ifquad xin Omega_{i+1},\\
0, qquad otherwise.\\
end{array}
ight. ag{4}$$
因此，我们需要获得的问题（1）的有限元解$u^n_h$可表示为
$$u^n_h=sumlimits_{j=1}^{M-1} u_j^n phi_j , ag{5}$$
其中$u^n_i$ 为问题（1）在点$(x_i,t_n)$处的数值解,将（5）式代入（3）式，我们可以获得求解问题（1）的有限元格式的代数方程组。.
$$Big( au A + B Big)u^n=Bu^{n-1}+ au widehat{f}. ag{6}$$
其中$u^n=[u_1^n,u_2^n,cdots,u_{M-1}^n]^T$, $A=(a_{ij})_{M-1 imes 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}$均为零，则通过计算，我们发现
$$egin{array}{rl}a_{i,i}&=int_Omega (phi_i)‘(phi_i)‘dx=sumlimits_{i=1}^{M-1} 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\\&=int_0^1 frac{1}{h}dxi+int_0^1 frac{1}{h} dxi\\&=frac{2}{h},qquad i=1,2,cdots, M-1. ag{7}end{array}$$
$$egin{array}{rl}a_{i-1,i}=a_{i,i-1}&=int_Omega (phi_i)‘(phi_{i-1})‘dx=sumlimits_{i=1}^{M-1} 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}dxi\\&=frac{-1}{h},qquad i=2,3,cdots,M-1. ag{8}end{array}$$
可知刚度矩阵$A$为一三对角矩阵，同理可知质量矩阵$B=(b_{ij})_{M-1 imes M-1}$也为三对角矩阵，其矩阵元素由如下计算获得
$$egin{array}{rl}b_{i,i}&=int_Omega (phi_i)^2dx=sumlimits_{i=1}^{M-1} 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^2dxi+hint_0^1 (1-xi)^2dxi\\&=frac{2h}{3},qquad i=1,2,cdots, M-1. ag{9}end{array}$$
$$egin{array}{rl}b_{i-1,i}=b_{i,i-1}&=int_Omega (phi_i)(phi_{i-1})dx=sumlimits_{i=1}^{M-1} int_{Omega_i} (phi_i)(phi_{i-1})dx\\&=int_{x_{i-1}}^{x_{i-1}} (phi_i)(phi_i)dx\\&=hint_0^1 xi(1-xi)dxi \\&=frac{h}{6},qquad i=2,3,cdots,M-1. ag{10}end{array}$$
到了这里，我们终于可以松一口气了，主要的东西我们已经准备好了，，，，等等，，，还有载荷向量$widehat{f}=[(f^n,phi_1),(f^n,phi_2),cdots,(f^n,phi_{M-1})]^T$还没计算，接下来给出载荷向量的计算
$$egin{array}Big( f^n,phi_i Big)&=int_Omega f^n,phi_i dx=sumlimits^{M-1}_{i=1} int_{Omega_{i}}f^n,phi_i dx\\&=int_{x_{i-1}}^{x_i} f(x)phi_idx+int_{x_i}^{x_{i+1}} f(x)phi_idx\\&=hint_0^1 f(x_{i-1}+hxi)xdx+hint_0^1 f(x_i+hxi)(1-xi)dxi. ag{11}end{array}$$
为计算以上的一重积分，我在程序中采用的是17点$Gauss$求积公式，绝对符合我们对算法精度的要求！不会$Gauss$积分公式的同学可以查阅文献[2]pp116页.具体程序实现如下
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
到这里，整个编程到此结束，下面我来展示程序运行结果
以上是关于快速学会有限元编程的主要内容，如果未能解决你的问题，请参考以下文章