[有限元方法阶段汇总篇] 有限元入门简单 1D 示例程序(Helmholtz 方程)
前言
一些链接
之前写过三篇基础的有限元基础入门级别的程序介绍,引起比较大的反响,很多人在 CSDN 私信和添加我微信讨论了相关问题,有限元入门见如下链接。
有限元方法入门:有限元方法简单的一维算例
有限元方法入门:有限元方法简单的二维算例(三角形剖分)
有限元方法入门:有限元方法简单的二维算例(矩形剖分)
同时,我也写了一些进阶的有限元介绍,也有非常多的人关注。进阶有限元介绍,见如下链接。
MATLAB 的 PDE 工具箱的简单使用
径向Kohn-Sham方程的谱有限元方法
曲面逼近和等参trace有限元
一维声子晶体的谱有限元方法
特征值问题的有限元MATLAB程序(一维
有限元方法基础入门教程(一维弹性问题的有限元程序)
单位分解有限元方法(PUFEM)
多尺度(有限元)降阶模型下的DNN方法
曲面偏微分方程:参数化有限元方法
神经网络和有限元方法
当然这部分链接并不是非常全,有其他需求的同学可以去我的博客里面再翻一翻。后续,我也会添加新的有限元相关文章,自然也不会在这里列出。但也建议先收藏本文。
动机
有很多人看了我有限元入门的三个链接,对于自己的问题,而不知道怎么使用我的程序。我可以负责地说,入门程序的框架非常完整,对于一些简单的问题,你要做的就是看懂我的程序,然后根据你自己的问题稍加修改即可。 不管是你想把程序改成 P2 元,或者是时间依赖的方程,比如说抛物方程,还是改成特征值问题,亦或是多尺度有限元,都是非常容易,因为这几个入门程序的框架还是很正规和传统的。譬如,对于相对应的抛物方程,你在时间上做一个隐式欧拉(也叫后向欧拉)差分,剩下的,就和求解椭圆方程,没有太大区别。
本文是对有限元方法三个入门算例的使用的一个补充说明。我接下来将以亥姆霍兹特征值问题为例,告诉大家如何将模板入门程序改写成能解决我们的问题的样子。
问题描述
关于特征值问题
我们知道,对于上述三个入门例子中,Lax Milgram 定理告诉我们,他的解是存在唯一的。而对于一般的特征值问题形式的方程,你再去验证 Lax Milgram 的条件,你会发现,不一定满足,这时候特征值问题对应的解就不一定是存在唯一的,可能不存在,存在也不一定唯一。比如,对于 − Δ u + b u = 0 -\Delta u+bu=0 −Δu+bu=0,只有 b b b 是非负的时候, Lax Milgram 定理才成立。负的拉普拉斯算子是一个半正定的算子,加上 0 本质边界条件,可以理解成了具有正定性。正定的算子,有正的特征值,所以当这里的 b b b取到拉普拉斯算子的特征值的时候,这个问题有非零解。这就是所谓的特征值问题。可见,Lax Milgram定理只是充分条件,而不是必要条件。
你也可以用代数方程的情况去类比推想方程的情况。对于 A x = b x Ax=bx Ax=bx, A A A 半正定。除零解外,若 b b b 是特征值,有一解多解的情况。若 b b b 是负的不是特征值,无解。考虑 b = 0 b=0 b=0,如果 0 0 0 是特征值,有线性解 b u + c bu+c bu+c,限定边界条件,这个线性解就定下来了。如果 0 0 0 不是特征值,无解。 0 0 0 特征值只能是单特征值,也就是说 b b b 取零,本质边界条件下,解都是唯一的。这段话比较拗口,可不看。
亥姆霍兹方程
亥姆霍茲方程 (英语: Helmholtz equation) 是一个描述电磁波的格圆偏微分方程, 以德国物理学家亥姆霍兹的名字命名。其基本形式如下:
(
∇
2
+
k
2
)
A
=
0
\left(\nabla^{2}+k^{2}\right) A=0
(∇2+k2)A=0
其中
∇
2
\nabla^{2}
∇2 是拉音拉斯算子,
k
k
k 是波数,
A
A
A 是振幅。
因为它和波动方程的关系,亥姆霍兹方程在物理学中电磁辐射、地震学和声学等相关研究领域里有着广泛应用。
程序与结果
问题一
求解问题:
d
2
u
d
x
2
+
k
u
(
x
)
=
l
(
x
)
\frac{d^{2} u}{d x^{2}}+k u(x)=l(x)
dx2d2u+ku(x)=l(x) 边界条件:
{
u
(
0
)
=
0
u
(
1
)
=
0
\left\{\begin{array}{l} u(0)=0 \\ u(1)=0 \end{array}\right.
{u(0)=0u(1)=0 定义
l
(
x
)
=
−
sin
(
π
x
)
l(x)=-\sin (\pi x)
l(x)=−sin(πx) 那么解析解为
u
(
x
)
=
1
π
2
−
k
sin
(
π
x
)
u(x)=\frac{1}{\pi^{2}-k} \sin (\pi x)
u(x)=π2−k1sin(πx) 当
k
=
0.5
π
2
k=0.5\pi^2
k=0.5π2时,解析解和数值解如图
当
k
=
π
2
k=\pi^2
k=π2 时,
当
k
=
1.0001
π
2
k=1.0001\pi^2
k=1.0001π2 时,
图表意思明确,今天圣诞,夜已深,废话就不多说了。对应的程序如下:
%% 求解一维 Helmholtz 方程的程序
% \laplace u + ku = l(x)
%% 一维有限元程序
clc; % 清空命令行窗口
clear; %清除工作空间
close all; % 关闭所有图像
for n = 2:2
%% 参数设置
%方程参数
k = 0.5*pi^2;
l = @(x) -sin(pi*x);
u = @(x) 1/(pi^2-k)*sin(pi*x);
L = 1; %定义单元长度,假定0为初始,L即为右边界
u_0 = 0; u_1 = 0; %定义第一类边界条件
% 数值参数
%n = 5;
numel = 2^n; %定义分割的单元数目
numnod = numel + 1; %节点个数比单元个数多1,节点编号等同于形函数编号
nel = 2; %每个单元的节点数目,即每个单元上有几个形函数参与作用,单元自由度
nsp = 2; %高斯点的个数(因为单元上的积分使用的是高斯积分公式),P1 元最多二阶精度,两个高斯点很够了
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,k);%计算刚度矩阵和质量矩阵
fe = elemforce(e,nel,nsp,coord,connect,l);
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:numnod
if (isempty(find(ebcdof == j, 1))) % 第j个点若不是固定点
fext(j) = fext(j) - bigk(j,n)*ebcval(i);%按固定点来处理右端项
end
end
%因为条件没用充分,刚度矩阵是不可逆的,我们要对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 = u(xsamp);
%% 绘图,可视化
plotflat = 1;
if plotflat==1
% if (numel > 20)
plot(coord,u_coeff,'-',xsamp,uexact,'--');
% else
% plot(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);
hold on;
drawnow;
end
%% 只计算 l_inf 误差
% u_ex = u(coord);
% error_Linf = max(abs(u_cal - u_ex))
error_Linf = -1;
u_cal_full = full(u_cal);
for x = linspace(0,1,numel*1000)
uh = interp1(coord,u_cal_full,x,'linear');
ut = u(x);
error_Linf = max(abs(uh-ut),error_Linf);
end
%format long;
vpa(error_Linf,3)
end
function 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;
else
error('GAUSS_WEIGHTS - Fatal error! Illegal value of n.')
end
return
end
function 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;
else
error('GAUSS_WEIGHTS - Fatal error! Illegal value of n.')
end
return
end
function [ke] = elemstiff(e,nel,nsp,coord,connect,k)
%输入单元编号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))];%计算两个不同的形函数在高斯点处的值,N'N为 4*4 矩阵,每个位置刚好表示 4 个被积函数在高斯点处的值
% xg = N*xe;%将高斯点映射到计算单元上去
A = N;% 表示两个基函数在高斯点处的值,因为两个基函数,所以有两个值,因为公式相同,直接引用N值
B = 1/Le*[-1 1];% 表示单元上两段线的斜率(导数值),即一般单元上的形函数导数值(于高斯点处)
ke = ke + weight(i)*(-B'*B+k*(A'*A))*detj;%计算单元刚度矩阵,导数值之间的作用直接在标准单元上积分就可以了,不用乘雅克比行列式
%等价于先在一般单元上求完导数之后,再映射到标准单元,这时候就要乘雅克比行列式了,因为两者的导数刚好差一个雅克比行列式
end
return
end
function [fe] = elemforce(e,nel,nsp,coord,connect,l)
%输入单元编号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:nsp
N = [ 0.5*(1 - xi(i)) 0.5*(1 + xi(i)) ];%容易验证 N.*[a b],表示的是高斯点 x(i) 到标准单元的映射
%N 也表示在标准单元上,两个单元形函数在高斯点上的值
xg = N*xe;%计算高斯点在一般单元上对应的位置,即将高斯点在标准单元映射会一般单元
bx = l(xg);%计算相应点的函数值
fe = fe + weight(i)*bx*N'*detj;%计算单元载荷,bx 是对应于高斯点的一个 f 值
end
return
end
问题二
试着求解以下的亥姆霍兹方程,
(
a
u
′
)
′
+
k
2
u
=
f
in
(
0
,
1
)
\left(a u^{\prime}\right)^{\prime}+k^{2} u=f \quad \text { in }(0,1)
(au′)′+k2u=f in (0,1)
边界条件是,
{
u
(
0
)
=
0
u
(
1
)
=
0
\left\{\begin{array}{l} u(0)=0 \\ u(1)=0 \end{array}\right.
{u(0)=0u(1)=0
这里的
a
(
x
)
:
=
1
+
1
2
sin
x
ε
a(x):=1+\frac{1}{2} \sin \frac{x}{\varepsilon}
a(x):=1+21sinεx
ε \varepsilon ε 是周期的, f f f 是光滑的。
测试你的代码对于
ε
=
1
2
n
,
n
=
0
,
1
,
2
,
3
,
4
\varepsilon=\frac{1}{2^{n}}, n=0,1,2,3,4
ε=2n1,n=0,1,2,3,4
通常地,为了得到精确解,你也需要改变你的网格尺寸:
h
∼
ε
/
10
h \sim \varepsilon / 10
h∼ε/10
所改进的代码如下:
%% 求解一维 Helmholtz 方程的程序
% (au')' + ku = l(x)
%% 一维有限元程序
clc; % 清空命令行窗口
clear; %清除工作空间
close all; % 关闭所有图像
for n = 0:4
%% 参数设置
%方程参数
k = 0.5*pi^2;
epsilon = 1/(2^n);
syms x;
a_syms = 1+0.5*sin(x/epsilon);
a = matlabFunction(a_syms);
u_syms = 1/(pi^2-k)*sin(pi*x);
u = matlabFunction(u_syms);
l_syms = diff(a_syms*diff(u_syms,x),x)+k*u_syms;
l = matlabFunction(l_syms);%右端项根据 u 来计算
L = 1; %定义单元长度,假定0为初始,L即为右边界
u_0 = 0; u_1 = 0; %定义第一类边界条件
% 数值参数
numel = round(10/epsilon); %定义分割的单元数目
numnod = numel + 1; %节点个数比单元个数多1,节点编号等同于形函数编号
nel = 2; %每个单元的节点数目,即每个单元上有几个形函数参与作用,单元自由度
nsp = 2; %高斯点的个数(因为单元上的积分使用的是高斯积分公式),P1 元最多二阶精度,两个高斯点很够了
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,k,a);%计算刚度矩阵和质量矩阵
fe = elemforce(e,nel,nsp,coord,connect,l);
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:numnod
if (isempty(find(ebcdof == j, 1))) % 第j个点若不是固定点
fext(j) = fext(j) - bigk(j,n)*ebcval(i);%按固定点来处理右端项
end
end
%因为条件没用充分,刚度矩阵是不可逆的,我们要对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 = u(xsamp);
%% 绘图,可视化
plotflat = 0;
if plotflat==1
% if (numel > 20)
plot(coord,u_coeff,'-',xsamp,uexact,'--');
% else
% plot(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);
hold on;
drawnow;
end
%% 只计算 l_inf 误差
% u_ex = u(coord);
% error_Linf = max(abs(u_cal - u_ex))
error_Linf = -1;
u_cal_full = full(u_cal);
for x = linspace(0,1,numel*1000)
uh = interp1(coord,u_cal_full,x,'linear');
ut = u(x);
error_Linf = max(abs(uh-ut),error_Linf);
end
%format long;
vpa(error_Linf,3)
end
function 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;
else
error('GAUSS_WEIGHTS - Fatal error! Illegal value of n.')
end
return
end
function 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;
else
error('GAUSS_WEIGHTS - Fatal error! Illegal value of n.')
end
return
end
function [fe] = elemforce(e,nel,nsp,coord,connect,l)
%输入单元编号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:nsp
N = [ 0.5*(1 - xi(i)) 0.5*(1 + xi(i)) ];%容易验证 N.*[a b],表示的是高斯点 x(i) 到标准单元的映射
%N 也表示在标准单元上,两个单元形函数在高斯点上的值
xg = N*xe;%计算高斯点在一般单元上对应的位置,即将高斯点在标准单元映射会一般单元
bx = l(xg);%计算相应点的函数值
fe = fe + weight(i)*bx*N'*detj;%计算单元载荷,bx 是对应于高斯点的一个 f 值
end
return
end
function [ke] = elemstiff(e,nel,nsp,coord,connect,k,a)
%输入单元编号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))];%计算两个不同的形函数在高斯点处的值,N'N为 4*4 矩阵,每个位置刚好表示 4 个被积函数在高斯点处的值
xg = N*xe;%将高斯点映射到计算单元上去
a_g = a(xg);%计算高斯点处的a值
A = N;% 表示两个基函数在高斯点处的值,因为两个基函数,所以有两个值,因为公式相同,直接引用N值
B = 1/Le*[-1 1];% 表示单元上两段线的斜率(导数值),即一般单元上的形函数导数值(于高斯点处)
ke = ke + weight(i)*(-a_g*(B'*B)+k*(A'*A))*detj;%计算单元刚度矩阵,导数值之间的作用直接在标准单元上积分就可以了,不用乘雅克比行列式
%等价于先在一般单元上求完导数之后,再映射到标准单元,这时候就要乘雅克比行列式了,因为两者的导数刚好差一个雅克比行列式
end
return
end
结果和前面的类似,我就不贴了,感兴趣的同学自己跑跑程序就知道了。
写在后面的话
怎么样?是不是很简单?主体框架和基本内容还是那三篇入门程序的,稍微改改就行了。这里的右端项,边值条件,变系数什么的,都是可以随便改的。
原来的时候,我回复大家的评论、私信和微信消息还比较勤快,后面渐渐地就回得少了甚至不回了。一方面我觉得我可能变得更忙了,另一方面,有些人加了我微信,我无偿回答完问题,就直接把我删了,这让我感觉心里有点拔凉拔凉的,所以,渐渐地就不想搭理一些人了。
再一个,如果你有需求,也可以让我帮忙写程序的,像这篇文章中这种的程序,我基本上收费也就在五百块钱,别的更复杂的,可以另议,不要讨价还价。问问题的话,你可以加我微信,先发个红包,再提出问题即可,这也从一定程度上能够预防那些过河拆桥的人的存在。
其实,如果你编程能力强的话,我更支持你们自力更生,丰衣足食。因为,你只要认真读懂了开头列出的几个入门程序,再做一个适配自己问题版本的有限元程序,是很快的。
由于一些令人失望的经验,我对于推广有限元方法的数学原理并不抱太乐观的态度。由于 Linux 操作系统和 C++ 这两个东西对于我们大多数同行和同学来说是两个巨大的困难,所以为了更照顾编程差的同学,我是用了大问题上运行速度极慢的 MATLAB 语言。如果您确实有研究上的需求,在弄懂原理的同时,需要对应地去动手编写程序,深耕半年,才有可能有所小成。在 当今这个日新月异的时代里,这个时间恐怕已经超出了大多数期望尽快发表 SCI 论文而能够付出的可接受时间尺度。
一般来说,在我的模板程序的帮助下,您写一个简单的偏微分方程的有限元程序的工作,在一到两个小时内完成是并不困难的。而且您的程序天生会具有相当大的可重用性:解一个新的问题常常只需要在旧程序上进行少许修改就能完成。它现在对于其能力所能够覆盖的问题,基本上能够比较完好的解决了,但是对于其能力不能覆盖的问题,是不可能在其现在的基本数据结构下进行简单的修改和补充就 能够解决的。
随着现在计算机硬件的计算能力的飞速发展,事实上来说,在过去的十几年间,我们的计算机的速度和容量提升了大约一千倍。我们作为来挖掘计算机的计算能力的研究人员,常常还在解决和十几年以前同样规模的问题写文章,因此事实上来说,我们的计算能力其实是退步了。我们依赖 Windows, 依赖 MATLAB,我们的研究生常常连两万阶的线性代数方程组都无法自己写程序解出来,我们越来越远离底层,无法理解计算数学到底在干什么和要干什么,成为漂浮在厚厚的一层虚无飘渺 之上的白痴。
期望阅读我的代码的过程,是让您从一个不同的视角来理解有限元作为一种计算方法的机会。这真是我 2020 年的最后一篇文章了,祝大家双旦快乐。
参考文献
FINITE ELEMENT ANALYSIS OF A ONE-DIMENSIONAL HELMHOLTZ EQUATION, JIANG, YIFAN DEPARTMENT OF ELECTRICAL ENGINEERING AND PHYSICS