有限元方法简单的一维算例

算例描述

我们对下述边值问题\label{eq1} {u′′(x)+u(x)=(1−π2)sinπx0≤x≤1u(0)=u(1)=0\left \{ \begin{aligned} & u''(x)+u(x) = (1-\pi^2)\text{sin}\pi x & 0\leq x \leq 1\\ & u(0) = u(1) = 0 & \end{aligned} \right.{u′′(x)+u(x)=(1−π2)sinπxu(0)=u(1)=00≤x≤1 考虑其变分问题，对其变分问题有限元离散并求解，并验证其收敛阶。
注：问题的真解为u(x)=sinπxu(x) = \text{sin}\pi xu(x)=sinπx。

变分问题

∀v∈H01\forall v \in H_0^1∀v∈H01，乘([eq1])式两边，使用格林公式，利用边界条件，易得\label{eq2}：
∫01u′v′+uvdx=∫01fvdx\int _0^1 u'v'+ uv dx = \int _0^1 fv dx∫01u′v′+uvdx=∫01fvdx 其中fff为方程中的右端项。令\label{eq3}

a(u,v)=∫01u′v′+uvdxf(v)=∫01fv\begin{aligned} &a(u,v) = \int _0^1 u'v'+ uv dx& \\ &f(v) = \int _0^1 fv& \end{aligned}a(u,v)=∫01u′v′+uvdxf(v)=∫01fv

容易证明，原问题([eq1])等价于变分问题: \label{eq4}{求u∈H01，使得a(u,v)=f(v)∀v∈H01\left \{ \begin{aligned} & \text{求}u \in H_0^1 \text{，使得}\\ & a(u,v) = f(v)&\forall v \in H_0^1 \end{aligned} \right.{求u∈H01，使得a(u,v)=f(v)∀v∈H01
事实上，在一定连续性的要求下，强解为弱解，弱解也是强解，二者等价。故求解问题([eq1])变为了求解问题([eq4])。更一般的变分问题，描述为：\label{eq5}
{求u∈V,s.t.a(u,v)=f(v)∀v∈V\left \{ \begin{aligned} & \text{求}u \in V ,s.t.&\\ & a(u,v) = f(v)&\forall v \in V& \end{aligned} \right.{求u∈V,s.t.a(u,v)=f(v)∀v∈V

有限元离散

问题转化

我们来考虑上述变分问题的有限维逼近，即构造VVV的有限维子空间Vh⊂VV_h \subset VVh⊂V，考虑如下的离散问题：\label{eq6}
{求uh∈Vh,s.t.a(uh,vh)=f(vh)∀vh∈Vh\left \{ \begin{aligned} & \text{求}u_h \in V_h ,s.t.&\\ & a(u_h,v_h) = f(v_h)&\forall v_h \in V_h& \end{aligned} \right.{求uh∈Vh,s.t.a(uh,vh)=f(vh)∀vh∈Vh

我们用问题([eq6])近似问题([eq5])，后者的解逼近前者的解。所以我们可以通过求([eq6])的解作为近似解。
设VhV_hVh空间中的一组基为{ϕi}i=1N\{\phi_i\}_{i=1}^N{ϕi}i=1N，若uh=Σuiϕiu_h = \Sigma u_i \phi_iuh=Σuiϕi（为了书写方便，不加说明，求和指标都为1到N），我们需要求的是组合系数uiu_iui，将uh=Σuiϕiu_h = \Sigma u_i \phi_iuh=Σuiϕi代入，并依次分别取vhv_hvh 为每个基函数，我们可以得到： \label{eq7}Σuia(ϕi,ϕj)=f(ϕj),j=1,2,...N.\Sigma u_i a(\phi_i,\phi_j)=f(\phi_j),j=1,2,...N.Σuia(ϕi,ϕj)=f(ϕj),j=1,2,...N.
这事实上就是一个线性方程组AU=FAU=FAU=F，其中A=(a(ϕj,ϕi)i,j=1N),U=(u1,...,uN)T,F=(f(ϕ1,...,ϕN))TA = (a(\phi_j,\phi_i)_{i,j=1}^N),U=(u_1,...,u_N)^T,F=(f(\phi_1,...,\phi_N))^TA=(a(ϕj,ϕi)i,j=1N),U=(u1,...,uN)T,F=(f(ϕ1,...,ϕN))T。那么，求解问题就转为了在VhV_hVh的约束下，求解这个线性方程组。

由我们选取vhv_hvh分别为ϕi\phi_iϕi可知，这个方程组的解对于原问题是必要而不充分的。但是在合适的条件之下，由Lax-Milgram定理可知，离散变分问题的解是存在唯一的。一般来说，我们在初边值条件的约束下，方程组的解存在唯一，那么这个唯一解必然是原问题的解。故而，对于一般的变分问题，我们通过离散化，最后可以通过求解([eq7])来近似。那么问题本质上就转化为了找一个合适的空间逼近，在这个空间中求基函数，和基函数之间的某种相互作用以及基函数和fff之间的作用，构建出方程组，最后求解方程组，得到逼近阶关于基函数的组合系数，最后得到原问题的一个逼近的数值解。

需要注意的是，由于Vh∈VV_h\in VVh∈V，我们对ϕi\phi_iϕi是有一定的要求的，或者说我们不管基函数是否属于VVV,我们最后对组合系数uiu_iui施以一定的要求（为满足边界条件），使得uh=Σuiϕi∈Vu_h = \Sigma u_i \phi_i \in Vuh=Σuiϕi∈V，具有相同的效力。

有限元三要素

我们现在来考虑单元上的插值问题。对于问题([eq1])，我们考虑其三要素(T,PT,ΣT)(T,P_T,\Sigma_T)(T,PT,ΣT)，其中，TTT为线段上的等距或者不等距分割，分割单元上的形函数空间为多项式集合PT=span{1,x}P_T = \text{span}\{1,x\}PT=span{1,x}，单元自由度ΣT={u(a1),u(a2)}\Sigma_T = \{u(a_1),u(a_2)\}ΣT={u(a1),u(a2)}
表示分割单元两端点的值。

适定性
一维线性插值，两点为零，插值函数恒为零，适定性满足。
求基函数
使用待定系数法，或者直接可以写出一维插值的节点基函数为：(x−a1)/(a1−a2)(x-a_1)/(a_1-a_2)(x−a1)/(a1−a2)和(x−a1)/(a2−a1)(x-a_1)/(a_2-a_1)(x−a1)/(a2−a1)。
连续性
显然，一维的时候节点处的值为u(ai)u(a_i)u(ai)，是连续的。

整体基函数和边值条件的处理

全局基函数

求得了在每个单元上的插值的节点基函数，对于每个节点，我们将其相关单元的相关基函数并成一个关于这个节点的全局基函数ϕi\phi_iϕi，有多少个节点，就有多少个基函数。
整理可得基函数如图([pic1])所示：

边值条件的处理

有了基函数，理论上就可以求解变分问题离散后转化成的方程组。但是在边值条件问题中，粗算出来的方程组的刚度矩阵KKK，就不是满秩的，那么解就不唯一。这是因为所求解的uh∈Vh⊂Vu_h \in V_h \subset Vuh∈Vh⊂V这个条件没满足，所以，一个思路是根据边界条件选取一个合适的基函数子集，另外，也可以最后根据边值条件来调整刚度矩阵KKK，我采取第二个思路。

调整的思路就是调整KU=FKU=FKU=F，将其等价转换，使得最后的解满足UUU的对应位置为固定的值，所谓对应位置，就是边界所对应的位置，所谓固定的值，就是边界值。简单地说，我们可以这样操作，直接令UUU的对应位置强制为边界值，然后将UUU中新的不知道的数视为未知数，构建新的方程组求解，本质上无非就是右端项减去边界条件乘以KKK中与其位置对应的列……假若，原来的KKK秩比起满秩少了m，那么补上m个边界条件，方程组的解应该是存在唯一的，认识到这一点，再加上一些线性代数的基本知识，事情就很明了了……就不再赘述……

数值实验和收敛阶

程序代码

基于以上的思想，编写了程序代码，直接粘贴如下，为了方便，直接将函数文件直接粘贴在主文件末尾了。程序的细节可以看注释，就不再详述。我基本上将每一句代码都打上了注释。

%% 一维有限元程序
clc;       % 清空命令行窗口
clear; %清除工作空间
close all; % 关闭所有图像
%% 参数设置
L = 1; %定义单元长度，假定0为初始，L即为右边界
numel = 5; %定义分割的单元数目
u_0 = 0; u_1 = 0; %定义第一类边界条件
numnod = numel + 1; %节点个数比单元个数多1，节点编号等同于形函数编号
nel = 2; %每个单元的节点数目,即每个单元上有几个形函数参与作用，单元自由度
nsp = 2; %高斯点的个数（因为单元上的积分使用的是高斯积分公式）
coord = linspace(0,L,numnod)'; %等分节点的坐标（为了方便，我这里采用等分的方式，事实上单元长度可以不一致），列向量
connect = [1:(numnod-1); 2:numnod]';%连接矩阵，表示每个单元周围的节点编号，也就是涉及的形函数编号，一行表示一个单元的情况，下标从1开始
ebcdof = [1,numnod];   % 强制性边界点的编号，本例子中是两端
ebcval = [u_0,u_1]; %强制性边界条件的取值，在第一个点的地方为u0,最后一个点为u_1
bigk = sparse(numnod,numnod); % 刚度矩阵[K]，初始化为0，使用稀疏矩阵存储
fext = sparse(numnod,1);      % 载荷向量{f},初始化为0%%  计算系数矩阵K和右端项f，单刚组装总刚
for e = 1:numel %按单元扫描ke = elemstiff(e,nel,nsp,coord,connect);fe = elemforce(e,nel,nsp,coord,connect);sctr = connect(e,:);bigk(sctr,sctr) = bigk(sctr,sctr) + ke;%按照位置组装总刚fext(sctr) = fext(sctr) + fe;
end
for i = 1:length(ebcdof)%边界条件的处理，处理总纲和载荷向量n = ebcdof(i);%强制性边界编号遍历for j = 1:numnodif (isempty(find(ebcdof == j, 1))) % 第j个点若不是固定点fext(j) = fext(j) - bigk(j,n)*ebcval(i);%按固定点来处理右端项endend%因为条件没用充分，刚度矩阵是不可逆的，我们要对K进行处理，即ui对应位置强制改成边界值%那么需要对方程组进行修正，即右端项要减去K的第n列乘un%置K的第n行第n列为零，nn位置为1，右端项第n位子为对应边界值%当然，我们也可以缩小K矩阵来处理，一个意思bigk(n,:) = 0.0;bigk(:,n) = 0.0;bigk(n,n) = 1.0;fext(n) = ebcval(i);
end
u_coeff = bigk\fext;%求出系数，事实上也是函数在对应点上的值,可用追赶法求，我这直接用内置函数了
u_cal = u_coeff;
%% 求精确解
nsamp = 101;
xsamp = linspace(0,L,nsamp);%100等分区间中间有100个数
uexact = exactsolution(xsamp);%% 绘图，可视化
if (numel > 20)plot(coord,u_coeff,'-',xsamp,uexact,'--');
elseplot(coord,u_coeff,'-',xsamp,uexact,'--',...coord,zeros(numnod,1),'mo','MarkerEdgeColor','k',...%k为黑色边界'MarkerFaceColor',[0.5 1 0.6],'MarkerSize',10);%圈圈大小为10%MarkerFaceColor：标记点内部的填充颜色 MarkerEdgeColor：标记点边缘的颜色，m紫红o圆圈
end
h = legend('FE solution','Exact solution','location','NorthEast');%摆放图例
set(h,'fontsize',9);%图例大小
title('Comparison of FE and Exact Solutions');%标题
xt = get(gca,'XTickLabel'); set(gca,'XTickLabel',xt,'fontsize',12);
%获取图像的横坐标tick，利用获取的tick，设置图像的横坐标标识
yt = get(gca,'XTickLabel'); set(gca,'XTickLabel',yt,'fontsize',12);
xlabel('$x$','Interpreter','latex','fontsize',16)%latex解释器下的横纵坐标名称显示
ylabel('$u^h , u$','Interpreter','latex','fontsize',16);%% 计算误差，我这里的误差只是使用节点处的值来估计，并没有计算单元积分作为准确的误差
%算收敛阶，大差不差了
u_ex = exactsolution(coord);
error_L1 = sum(abs(u_cal - u_ex))/(numnod-1) %算每个节点的平均误差
error_L2 = sqrt(sum((u_cal - u_ex).^2)/(numnod-1))
error_Linf = max(abs(u_cal - u_ex))function [ke] = elemstiff(e,nel,nsp,coord,connect)
%输入单元编号e，单元自由度nel，高斯点数目nsp，等分节点的坐标coord和连接矩阵connect
%输出单元刚度矩阵，是叫单元刚度矩阵吗？
ke    = zeros(nel,nel); %单刚初始化
nodes = connect(e,:);%单元相关形函数（节点）编号
xe    = coord(nodes);%相关节点的坐标
Le    = xe(nel) - xe(1);%单元长度，表示一种细度
detj  = Le/2;%雅克比行列式（一维）即为长度和标准单元长度的比值
xi = gauss_points(nsp);%选取高斯积分点【-1 1】上的
weight = gauss_weights(nsp);%高斯积分点对应的权重
for i = 1:nsp%按高斯点来进行积分计算，不同形函数之间的相互作用用列向量乘行向量来体现N = [ 0.5*(1 - xi(i))  0.5*(1 + xi(i))];%从标准单元映射到一般单元对于每个节点所需要的权重% xg = N*xe;%将高斯点映射到计算单元上去A = N;% 表示两个基函数在高斯点处的值，因为两个基函数，所以有两个值，因为公式相同，直接引用N值B = 1/Le*[-1 1];% 表示单元上两段线的斜率（导数值），即一般单元上的形函数导数值（于高斯点处）ke = ke + weight(i)*(-B'*B+A'*A)*detj;%计算单元刚度矩阵
endreturn
endfunction [fe] = elemforce(e,nel,nsp,coord,connect)
%输入单元编号3，单元自由度nel，高斯点数目nsp，等分坐标点coord和连接矩阵connect
%输入单元载荷向量
fe = zeros(nel,1);%初始化单元载荷向量为0
nodes = connect(e,:);%单元相关节点的编号
xe = coord(nodes); % 单元相关节点坐标
Le  = xe(nel) - xe(1);%单元长度
detj = Le/2;%雅克比行列式
xi = gauss_points(nsp);%高斯点
weight = gauss_weights(nsp);%高斯点对应的高斯权重
for i = 1:nspN = [ 0.5*(1 - xi(i))  0.5*(1 + xi(i)) ];%从标准单元映射到一般单元对于每个节点所需要的权重xg = N*xe;%计算高斯点在一般单元上对应的位置，即将高斯点在标准单元映射会一般单元bx = fun(xg);%计算相应点的函数值fe = fe + weight(i)*bx*N'*detj;%计算单元载荷
endreturn
end
function bx = fun(x)
bx = (1-pi^2)*sin(pi*x);
endfunction x = guass_points(n)
%
%  function x = gauss_points(n)
%
%  For 1 <= n <= 20, returns the abscissas x of an n point
%  Gauss-Legendre quadrature rule over the interval [-1,1].
%x = ones(1,n);if ( n == 1 )x(1) = 0.0;elseif ( n == 2 )x(1) = - 0.577350269189625764509148780502;x(2) =   0.577350269189625764509148780502;elseif ( n == 3 )x(1) = - 0.774596669241483377035853079956;x(2) =  0.000000000000000;x(3) =   0.774596669241483377035853079956;elseif ( n == 4 )x(1) = - 0.861136311594052575223946488893;x(2) = - 0.339981043584856264802665759103;x(3) =   0.339981043584856264802665759103;x(4) =   0.861136311594052575223946488893;elseif ( n == 5 )x(1) = - 0.906179845938663992797626878299;x(2) = - 0.538469310105683091036314420700;x(3) =   0.0;x(4) =   0.538469310105683091036314420700;x(5) =   0.906179845938663992797626878299;elseif ( n == 6 )x(1) = - 0.932469514203152027812301554494;x(2) = - 0.661209386466264513661399595020;x(3) = - 0.238619186083196908630501721681;x(4) =   0.238619186083196908630501721681;x(5) =   0.661209386466264513661399595020;x(6) =   0.932469514203152027812301554494;elseif ( n == 7 )x(1) = - 0.949107912342758524526189684048;x(2) = - 0.741531185599394439863864773281;x(3) = - 0.405845151377397166906606412077;x(4) =   0.0;x(5) =   0.405845151377397166906606412077;x(6) =   0.741531185599394439863864773281;x(7) =   0.949107912342758524526189684048;elseif ( n == 8 )x(1) = - 0.960289856497536231683560868569;x(2) = - 0.796666477413626739591553936476;x(3) = - 0.525532409916328985817739049189;x(4) = - 0.183434642495649804939476142360;x(5) =   0.183434642495649804939476142360;x(6) =   0.525532409916328985817739049189;x(7) =   0.796666477413626739591553936476;x(8) =   0.960289856497536231683560868569;elseif ( n == 9 )x(1) = - 0.968160239507626089835576202904;x(2) = - 0.836031107326635794299429788070;x(3) = - 0.613371432700590397308702039341;x(4) = - 0.324253423403808929038538014643;x(5) =   0.0;x(6) =   0.324253423403808929038538014643;x(7) =   0.613371432700590397308702039341;x(8) =   0.836031107326635794299429788070;x(9) =   0.968160239507626089835576202904;elseif ( n == 10 )x(1) =  - 0.973906528517171720077964012084;x(2) =  - 0.865063366688984510732096688423;x(3) =  - 0.679409568299024406234327365115;x(4) =  - 0.433395394129247290799265943166;x(5) =  - 0.148874338981631210884826001130;x(6) =    0.148874338981631210884826001130;x(7) =    0.433395394129247290799265943166;x(8) =    0.679409568299024406234327365115;x(9) =    0.865063366688984510732096688423;x(10) =   0.973906528517171720077964012084;elseif (n == 11)x(1)=-0.978228658146056992803938001123;x(2)=-0.887062599768095299075157769304;x(3)=-0.730152005574049324093416252031;x(4)=-0.519096129206811815925725669459;x(5)=-0.269543155952344972331531985401;x(6)= 0;x(7)= 0.269543155952344972331531985401;x(8)= 0.519096129206811815925725669459;x(9)= 0.730152005574049324093416252031;x(10)= 0.887062599768095299075157769304;x(11)= 0.978228658146056992803938001123;elseif (n == 12)x(1)=-0.981560634246719250690549090149;x(2)=-0.904117256370474856678465866119;x(3)=-0.769902674194304687036893833213;x(4)=-0.587317954286617447296702418941;x(5)=-0.367831498998180193752691536644;x(6)=-0.125233408511468915472441369464;x(7)= 0.125233408511468915472441369464;x(8)= 0.367831498998180193752691536644;x(9)= 0.587317954286617447296702418941;x(10)= 0.769902674194304687036893833213;x(11)= 0.904117256370474856678465866119;x(12)= 0.981560634246719250690549090149;elseif (n == 13)x(1)=-0.984183054718588149472829448807;x(2)=-0.917598399222977965206547836501;x(3)=-0.801578090733309912794206489583;x(4)=-0.642349339440340220643984606996;x(5)=-0.448492751036446852877912852128;x(6)=-0.230458315955134794065528121098;x(7)= 0;x(8)= 0.230458315955134794065528121098;x(9)= 0.448492751036446852877912852128;x(10)= 0.642349339440340220643984606996;x(11)= 0.801578090733309912794206489583;x(12)= 0.917598399222977965206547836501;x(13)= 0.984183054718588149472829448807;elseif (n == 14)x(1)=-0.986283808696812338841597266704;x(2)=-0.928434883663573517336391139378;x(3)=-0.827201315069764993189794742650;x(4)=-0.687292904811685470148019803019;x(5)=-0.515248636358154091965290718551;x(6)=-0.319112368927889760435671824168;x(7)=-0.108054948707343662066244650220;x(8)= 0.108054948707343662066244650220;x(9)= 0.319112368927889760435671824168;x(10)= 0.515248636358154091965290718551;x(11)= 0.687292904811685470148019803019;x(12)= 0.827201315069764993189794742650;x(13)= 0.928434883663573517336391139378;x(14)= 0.986283808696812338841597266704;elseif (n == 15)x(1)=-0.987992518020485428489565718587;x(2)=-0.937273392400705904307758947710;x(3)=-0.848206583410427216200648320774;x(4)=-0.724417731360170047416186054614;x(5)=-0.570972172608538847537226737254;x(6)=-0.394151347077563369897207370981;x(7)=-0.201194093997434522300628303395;x(8)= 0;x(9)= 0.201194093997434522300628303395;x(10)= 0.394151347077563369897207370981;x(11)= 0.570972172608538847537226737254;x(12)= 0.724417731360170047416186054614;x(13)= 0.848206583410427216200648320774;x(14)= 0.937273392400705904307758947710;x(15)= 0.987992518020485428489565718587;elseif (n == 16)x(1)=-0.989400934991649932596154173450;x(2)=-0.944575023073232576077988415535;x(3)=-0.865631202387831743880467897712;x(4)=-0.755404408355003033895101194847;x(5)=-0.617876244402643748446671764049;x(6)=-0.458016777657227386342419442984;x(7)=-0.281603550779258913230460501460;x(8)=-0.0950125098376374401853193354250;x(9)= 0.0950125098376374401853193354250;x(10)= 0.281603550779258913230460501460;x(11)= 0.458016777657227386342419442984;x(12)= 0.617876244402643748446671764049;x(13)= 0.755404408355003033895101194847;x(14)= 0.865631202387831743880467897712;x(15)= 0.944575023073232576077988415535;x(16)= 0.989400934991649932596154173450;elseif (n == 17)x(1)=-0.990575475314417335675434019941;x(2)=-0.950675521768767761222716957896;x(3)=-0.880239153726985902122955694488;x(4)=-0.781514003896801406925230055520;x(5)=-0.657671159216690765850302216643;x(6)=-0.512690537086476967886246568630;x(7)=-0.351231763453876315297185517095;x(8)=-0.178484181495847855850677493654;x(9)= 0;x(10)= 0.178484181495847855850677493654;x(11)= 0.351231763453876315297185517095;x(12)= 0.512690537086476967886246568630;x(13)= 0.657671159216690765850302216643;x(14)= 0.781514003896801406925230055520;x(15)= 0.880239153726985902122955694488;x(16)= 0.950675521768767761222716957896;x(17)= 0.990575475314417335675434019941;elseif (n == 18)x(1)=-0.991565168420930946730016004706;x(2)=-0.955823949571397755181195892930;x(3)=-0.892602466497555739206060591127;x(4)=-0.803704958972523115682417455015;x(5)=-0.691687043060353207874891081289;x(6)=-0.559770831073947534607871548525;x(7)=-0.411751161462842646035931793833;x(8)=-0.251886225691505509588972854878;x(9)=-0.0847750130417353012422618529358;x(10)= 0.0847750130417353012422618529358;x(11)= 0.251886225691505509588972854878;x(12)= 0.411751161462842646035931793833;x(13)= 0.559770831073947534607871548525;x(14)= 0.691687043060353207874891081289;x(15)= 0.803704958972523115682417455015;x(16)= 0.892602466497555739206060591127;x(17)= 0.955823949571397755181195892930;x(18)= 0.991565168420930946730016004706;elseif (n == 19)x(1)=-0.992406843843584403189017670253;x(2)=-0.960208152134830030852778840688;x(3)=-0.903155903614817901642660928532;x(4)=-0.822714656537142824978922486713;x(5)=-0.720966177335229378617095860824;x(6)=-0.600545304661681023469638164946;x(7)=-0.464570741375960945717267148104;x(8)=-0.316564099963629831990117328850;x(9)=-0.160358645640225375868096115741;x(10)= 0;x(11)= 0.160358645640225375868096115741;x(12)= 0.316564099963629831990117328850;x(13)= 0.464570741375960945717267148104;x(14)= 0.600545304661681023469638164946;x(15)= 0.720966177335229378617095860824;x(16)= 0.822714656537142824978922486713;x(17)= 0.903155903614817901642660928532;x(18)= 0.960208152134830030852778840688;x(19)= 0.992406843843584403189017670253;elseif (n == 20)x(1)=-0.993128599185094924786122388471;x(2)=-0.963971927277913791267666131197;x(3)=-0.912234428251325905867752441203;x(4)=-0.839116971822218823394529061702;x(5)=-0.746331906460150792614305070356;x(6)=-0.636053680726515025452836696226;x(7)=-0.510867001950827098004364050955;x(8)=-0.373706088715419560672548177025;x(9)=-0.227785851141645078080496195369;x(10)=-0.0765265211334973337546404093988;x(11)= 0.0765265211334973337546404093988;x(12)= 0.227785851141645078080496195369;x(13)= 0.373706088715419560672548177025;x(14)= 0.510867001950827098004364050955;x(15)= 0.636053680726515025452836696226;x(16)= 0.746331906460150792614305070356;x(17)= 0.839116971822218823394529061702;x(18)= 0.912234428251325905867752441203;x(19)= 0.963971927277913791267666131197;x(20)= 0.993128599185094924786122388471;elseerror('GAUSS_WEIGHTS - Fatal error! Illegal value of n.')endreturn
endfunction w = gauss_weights(n);
%
%  function w = gauss_weights(n);
%
%  For 1 <= n <= 20, returns the weights W of an
%  n point Gauss-Legendre quadrature rule over the interval [-1,1].
%w = ones(1,n);if ( n == 1 )w(1) = 2.0;elseif ( n == 2 )w(1) =  1.0;w(2) =  w(1);elseif ( n == 3 )w(1) =  0.5555555555555555555555555555565;w(2) =  0.8888888888888888888888888888889;w(3) =  0.5555555555555555555555555555565;elseif ( n == 4 )w(1) = 0.347854845137453857373063949222;w(2) = 0.652145154862546142626936050778;w(3) = 0.652145154862546142626936050778;w(4) = 0.347854845137453857373063949222;elseif ( n == 5 )w(1) = 0.236926885056189087514264040720;w(2) = 0.478628670499366468041291514836;w(3) = 0.568888888888888888888888888889;w(4) = 0.478628670499366468041291514836;w(5) = 0.236926885056189087514264040720;elseif ( n == 6 )w(1) = 0.171324492379170345040296142173;w(2) = 0.360761573048138607569833513838;w(3) = 0.467913934572691047389870343990;w(4) = 0.467913934572691047389870343990;w(5) = 0.360761573048138607569833513838;w(6) = 0.171324492379170345040296142173;elseif ( n == 7 )w(1) = 0.129484966168869693270611432679;w(2) = 0.279705391489276667901467771424;w(3) = 0.381830050505118944950369775489;w(4) = 0.417959183673469387755102040816;w(5) = 0.381830050505118944950369775489;w(6) = 0.279705391489276667901467771424;w(7) = 0.129484966168869693270611432679;elseif ( n == 8 )w(1) = 0.101228536290376259152531354310;w(2) = 0.222381034453374470544355994426;w(3) = 0.313706645877887287337962201987;w(4) = 0.362683783378361982965150449277;w(5) = 0.362683783378361982965150449277;w(6) = 0.313706645877887287337962201987;w(7) = 0.222381034453374470544355994426;w(8) = 0.101228536290376259152531354310;elseif ( n == 9 )w(1) = 0.0812743883615744119718921581105;w(2) = 0.180648160694857404058472031243;w(3) = 0.260610696402935462318742869419;w(4) = 0.312347077040002840068630406584;w(5) = 0.330239355001259763164525069287;w(6) = 0.312347077040002840068630406584;w(7) = 0.260610696402935462318742869419;w(8) = 0.180648160694857404058472031243;w(9) = 0.0812743883615744119718921581105;elseif ( n == 10 )w(1) =  0.0666713443086881375935688098933;w(2) =  0.149451349150580593145776339658;w(3) =  0.219086362515982043995534934228;w(4) =  0.269266719309996355091226921569;w(5) =  0.295524224714752870173892994651;w(6) =  0.295524224714752870173892994651;w(7) =  0.269266719309996355091226921569;w(8) =  0.219086362515982043995534934228;w(9) =  0.149451349150580593145776339658;w(10) = 0.0666713443086881375935688098933;elseif ( n == 11 )w(1)= 0.0556685671161736664827537204425;w(2)= 0.125580369464904624634694299224;w(3)= 0.186290210927734251426097641432;w(4)= 0.233193764591990479918523704843;w(5)= 0.262804544510246662180688869891;w(6)= 0.272925086777900630714483528336;w(7)= 0.262804544510246662180688869891;w(8)= 0.233193764591990479918523704843;w(9)= 0.186290210927734251426097641432;w(10)= 0.125580369464904624634694299224;w(11)= 0.0556685671161736664827537204425;elseif ( n == 12 )w(1)= 0.0471753363865118271946159614850;w(2)= 0.106939325995318430960254718194;w(3)= 0.160078328543346226334652529543;w(4)= 0.203167426723065921749064455810;w(5)= 0.233492536538354808760849898925;w(6)= 0.249147045813402785000562436043;w(7)= 0.249147045813402785000562436043;w(8)= 0.233492536538354808760849898925;w(9)= 0.203167426723065921749064455810;w(10)= 0.160078328543346226334652529543;w(11)= 0.106939325995318430960254718194;w(12)= 0.0471753363865118271946159614850;elseif ( n == 13 )w(1)= 0.0404840047653158795200215922010;w(2)= 0.0921214998377284479144217759538;w(3)= 0.138873510219787238463601776869;w(4)= 0.178145980761945738280046691996;w(5)= 0.207816047536888502312523219306;w(6)= 0.226283180262897238412090186040;w(7)= 0.232551553230873910194589515269;w(8)= 0.226283180262897238412090186040;w(9)= 0.207816047536888502312523219306;w(10)= 0.178145980761945738280046691996;w(11)= 0.138873510219787238463601776869;w(12)= 0.0921214998377284479144217759538;w(13)= 0.0404840047653158795200215922010;elseif ( n == 14 )w(1)= 0.0351194603317518630318328761382;w(2)= 0.0801580871597602098056332770629;w(3)= 0.121518570687903184689414809072;w(4)= 0.157203167158193534569601938624;w(5)= 0.185538397477937813741716590125;w(6)= 0.205198463721295603965924065661;w(7)= 0.215263853463157790195876443316;w(8)= 0.215263853463157790195876443316;w(9)= 0.205198463721295603965924065661;w(10)= 0.185538397477937813741716590125;w(11)= 0.157203167158193534569601938624;w(12)= 0.121518570687903184689414809072;w(13)= 0.0801580871597602098056332770629;w(14)= 0.0351194603317518630318328761382;elseif ( n == 15 )w(1)= 0.0307532419961172683546283935772;w(2)= 0.0703660474881081247092674164507;w(3)= 0.107159220467171935011869546686;w(4)= 0.139570677926154314447804794511;w(5)= 0.166269205816993933553200860481;w(6)= 0.186161000015562211026800561866;w(7)= 0.198431485327111576456118326444;w(8)= 0.202578241925561272880620199968;w(9)= 0.198431485327111576456118326444;w(10)= 0.186161000015562211026800561866;w(11)= 0.166269205816993933553200860481;w(12)= 0.139570677926154314447804794511;w(13)= 0.107159220467171935011869546686;w(14)= 0.0703660474881081247092674164507;w(15)= 0.0307532419961172683546283935772;elseif ( n == 16 )w(1)= 0.0271524594117540948517805724560;w(2)= 0.0622535239386478928628438369944;w(3)= 0.0951585116824927848099251076022;w(4)= 0.124628971255533872052476282192;w(5)= 0.149595988816576732081501730547;w(6)= 0.169156519395002538189312079030;w(7)= 0.182603415044923588866763667969;w(8)= 0.189450610455068496285396723208;w(9)= 0.189450610455068496285396723208;w(10)= 0.182603415044923588866763667969;w(11)= 0.169156519395002538189312079030;w(12)= 0.149595988816576732081501730547;w(13)= 0.124628971255533872052476282192;w(14)= 0.0951585116824927848099251076022;w(15)= 0.0622535239386478928628438369944;w(16)= 0.0271524594117540948517805724560;elseif ( n == 17 )w(1)= 0.0241483028685479319601100262876;w(2)= 0.0554595293739872011294401653582;w(3)= 0.0850361483171791808835353701911;w(4)= 0.111883847193403971094788385626;w(5)= 0.135136368468525473286319981702;w(6)= 0.154045761076810288081431594802;w(7)= 0.168004102156450044509970663788;w(8)= 0.176562705366992646325270990113;w(9)= 0.179446470356206525458265644262;w(10)= 0.176562705366992646325270990113;w(11)= 0.168004102156450044509970663788;w(12)= 0.154045761076810288081431594802;w(13)= 0.135136368468525473286319981702;w(14)= 0.111883847193403971094788385626;w(15)= 0.0850361483171791808835353701911;w(16)= 0.0554595293739872011294401653582;w(17)= 0.0241483028685479319601100262876;elseif ( n == 18 )w(1)= 0.0216160135264833103133427102665;w(2)= 0.0497145488949697964533349462026;w(3)= 0.0764257302548890565291296776166;w(4)= 0.100942044106287165562813984925;w(5)= 0.122555206711478460184519126800;w(6)= 0.140642914670650651204731303752;w(7)= 0.154684675126265244925418003836;w(8)= 0.164276483745832722986053776466;w(9)= 0.169142382963143591840656470135;w(10)= 0.169142382963143591840656470135;w(11)= 0.164276483745832722986053776466;w(12)= 0.154684675126265244925418003836;w(13)= 0.140642914670650651204731303752;w(14)= 0.122555206711478460184519126800;w(15)= 0.100942044106287165562813984925;w(16)= 0.0764257302548890565291296776166;w(17)= 0.0497145488949697964533349462026;w(18)= 0.0216160135264833103133427102665;elseif ( n == 19 )w(1)= 0.0194617882297264770363120414644;w(2)= 0.0448142267656996003328381574020;w(3)= 0.0690445427376412265807082580060;w(4)= 0.0914900216224499994644620941238;w(5)= 0.111566645547333994716023901682;w(6)= 0.128753962539336227675515784857;w(7)= 0.142606702173606611775746109442;w(8)= 0.152766042065859666778855400898;w(9)= 0.158968843393954347649956439465;w(10)= 0.161054449848783695979163625321;w(11)= 0.158968843393954347649956439465;w(12)= 0.152766042065859666778855400898;w(13)= 0.142606702173606611775746109442;w(14)= 0.128753962539336227675515784857;w(15)= 0.111566645547333994716023901682;w(16)= 0.0914900216224499994644620941238;w(17)= 0.0690445427376412265807082580060;w(18)= 0.0448142267656996003328381574020;w(19)= 0.0194617882297264770363120414644;elseif ( n == 20 )w(1)= 0.0176140071391521183118619623519;w(2)= 0.0406014298003869413310399522749;w(3)= 0.0626720483341090635695065351870;w(4)= 0.0832767415767047487247581432220;w(5)= 0.101930119817240435036750135480;w(6)= 0.118194531961518417312377377711;w(7)= 0.131688638449176626898494499748;w(8)= 0.142096109318382051329298325067;w(9)= 0.149172986472603746787828737002;w(10)= 0.152753387130725850698084331955;w(11)= 0.152753387130725850698084331955;w(12)= 0.149172986472603746787828737002;w(13)= 0.142096109318382051329298325067;w(14)= 0.131688638449176626898494499748;w(15)= 0.118194531961518417312377377711;w(16)= 0.101930119817240435036750135480;w(17)= 0.0832767415767047487247581432220;w(18)= 0.0626720483341090635695065351870;w(19)= 0.0406014298003869413310399522749;w(20)= 0.0176140071391521183118619623519;elseerror('GAUSS_WEIGHTS - Fatal error! Illegal value of n.')endreturn
endfunction uexact = exactsolution(xg)
uexact = sin(pi*xg);return
end

这里对于积分的计算(a(ϕi,ϕj)a(\phi_i,\phi_j)a(ϕi,ϕj)和f(ϕj)f(\phi_j)f(ϕj))我是采用扫描单元的方式，即在每个单元上计算基函数之间的相互作用并分配到相对应的全局基函数的相互作用上，计算每个单元上的基函数和右端项的作用，分配到与之相关的全局基函数和右端项的作用上，扫描全部单元，进行累加。

值得注意的是，我这里求积分没有利用基函数的特殊性，而是将其考虑为一般的基函数，提高了程序的泛化能力。另一方面，我求积分，没有利用matlab的内置函数，而是采用高斯积分公式来求积分，主要考虑到精度问题和速度问题。

程序中，我使用了很大篇幅去吧高斯点列出来，而不使用程序即用即算，主要是为了节省时间，以代码的直白换取时间的提高。这也是程序的一大亮点。事实上高斯点可以用如下代码求得，不过是最高阶正交多项式的零点。

p(1)=sym(a);p(2)=sym(b);
for i=3:k+1
syms x;
delta(i)=int(w_x*x*p(i-1)^2,x,a,b)/int(w_x*p(i-1)^2,x,a,b);
gamma2(i)=int(w_x*p(i-1)^2,x,a,b)/int(w_x*p(i-2)^2,x,a,b);
gamma2(3)=0;
p(i)=(x-delta(i))*p(i-1)-gamma2(i)*p(i-2);
end
p(1:k)=p(2:k+1);
p(k+1)=[];
p%正交多项式集合
%p=expand(p)
%做w(x)=1的k个正交多项式存放到数组p中；
%%
coe=sym2poly(p(k));
xn=sort(roots(coe)');%最高阶正交多项式零点

美中不足的是，代码中求解方程组部分，我直接使用matlab的右除算子\\backslash\，对于这种三对角的稀疏矩阵，其实可以考虑使用追赶法等。

误差和收敛阶分析

高斯求积公式的代数精确度是2n-1，所以对于一次多项式，取高斯点个数为1就很够了，也就是说只要看中点的值，就能得到一次多项式准确的积分值。但是考虑到舍入误差，我们不妨定义高斯点个数为2。定义单元长度为1，分割的单元数目为k，结果如下图所示。收敛阶的计算是基于L2L_2L2误差的，使用其它度量结果差不多。

由此可见，一维线性元的收敛阶是2。刚开始的时候比较小，只有1.4870，后面会趋于稳定，稳定在2附近。为什么会是2阶收敛呢？搞不懂的可以看看误差分析这一块。

有限元方法入门：有限元方法简单的一维算例相关推荐

方法入门_方法的调用
方法的调用方法在定义完毕后,方法不会自己运行,必须被调用才能执行,我们可以在主方法main中来调用我们自己定义好的方法.在主方法中,直接写要调用的方法名字就可以调用了. public static ...
有限体积法（13）——SIMPLE算法的一维算例
算例模型本文将使用SIMPLE算法来计算一个简单算例,以展示SIMPLE算法的具体事实过程.计算模型是一个平面喷嘴结构,尺寸数据如下图,假设流动是稳态无粘不可压的,计算速度和压力分布. 参数: 流体 ...
方法入门_方法的定义
方法的定义定义格式: 修饰符返回值类型方法名 (参数列表){代码...return ; } 定义格式解释: 修饰符: 目前固定写法 public static . 返回值类型: 目前固定写法 v ...
综合能源系统通用建模及规划方法研究-算例-笔记
综合能源系统通用建模及规划方法研究--算例分析区域综合能源系统规划对象包括能量枢纽和区域能源配网,前者的规划需综合考虑系统结构和设备型号,为区域综合能源系统规划的核心内容,为了更直接说明本文提出 ...
有限元方法入门：有限元方法简单的二维算例（三角形剖分）
有限元方法简单的二维算例(三角形剖分) 文章目录有限元方法简单的二维算例(三角形剖分) 算例描述变分问题有限元离散问题转化有限元三要素参考单元与一般单元一般单元上的形函数一般单元上的积 ...
有限元方法入门：有限元方法简单的二维算例（矩形剖分）
#有限元方法简单的二维算例(矩形剖分) 算例描述我们对下述椭圆边值问题 \label{eq1}{−Δu=fu∣∂Ω=0\left \{ \begin{aligned} & -\Delta u ...
二维三角元有限元方法matlab,有限元C++编程实践.doc
有限元C++编程实践.doc 基于有限元算法的编程实践学号:2011043010031 姓名:廖校毅电子科技大学物理电子学院[摘要]本文就有限元算法在电磁场分析中的应用展开研究与实践,从电磁场的 M ...
有限元方法基础-以二维拉普拉斯方程为例(附程序)
文章目录前言问题描述问题区域偏微分方程变分问题(弱形式) 有限元离散二维线性有限元 : 参考基函数 2D linear finite element : affine mapping 有限 ...
python划分有限元网格_有限元网格划分的基本原则及通用方法（有限元科技内参）...
的划分一方面要考虑对各物体几何形状的准确描述,另一方面也要考虑变形梯度的准确描述.为正确.合理地建立有限元模型,这里介绍划分网格时应考虑的一些基本原则. (1)网格数量网格数量直接影响计算精度和计算 ...
tensorflow打印模型结构_钢结构模型3D打印与有限元网格的融合方法
作者:魏鲁双刘尚蔚王颖魏群华北水利水电大学钢结构与工程研究院中国科学院大学人工智能学院摘要薄壁结构的3D打印STL文件是单侧外表面三角面网络,而钢结构模型体系的3D打印数据图形是由 ...

有限元方法入门：有限元方法简单的一维算例