[有限元方法階段彙總篇] 有限元入門簡單 1D 示例程式（Helmholtz 方程）

前言

一些連結

之前寫過三篇基礎的有限元基礎入門級別的程式介紹，引起比較大的反響，很多人在 CSDN 私信和新增我微信討論了相關問題，有限元入門見如下連結。

有限元方法入門：有限元方法簡單的一維算例
有限元方法入門：有限元方法簡單的二維算例（三角形剖分）
有限元方法入門：有限元方法簡單的二維算例（矩形剖分）

同時，我也寫了一些進階的有限元介紹，也有非常多的人關注。進階有限元介紹，見如下連結。

MATLAB 的 PDE 工具箱的簡單使用
徑向Kohn-Sham方程的譜有限元方法
曲面逼近和等參trace有限元
一維聲子晶體的譜有限元方法
特徵值問題的有限元MATLAB程式（一維
有限元方法基礎入門教程（一維彈性問題的有限元程式）
單位分解有限元方法(PUFEM)
多尺度（有限元）降階模型下的DNN方法
曲面偏微分方程：引數化有限元方法
神經網路和有限元方法

當然這部分連結並不是非常全，有其他需求的同學可以去我的部落格裡面再翻一翻。後續，我也會新增新的有限元相關文章，自然也不會在這裡列出。但也建議先收藏本文。

動機

有很多人看了我有限元入門的三個連結，對於自己的問題，而不知道怎麼使用我的程式。我可以負責地說，入門程式的框架非常完整，對於一些簡單的問題，你要做的就是看懂我的程式，然後根據你自己的問題稍加修改即可。 不管是你想把程式改成 P2 元，或者是時間依賴的方程，比如說拋物方程，還是改成特徵值問題，亦或是多尺度有限元，都是非常容易，因為這幾個入門程式的框架還是很正規和傳統的。譬如，對於相對應的拋物方程，你在時間上做一個隱式尤拉（也叫後向尤拉）差分，剩下的，就和求解橢圓方程，沒有太大區別。

本文是對有限元方法三個入門算例的使用的一個補充說明。我接下來將以亥姆霍茲特徵值問題為例，告訴大家如何將模板入門程式改寫成能解決我們的問題的樣子。

問題描述

關於特徵值問題

我們知道，對於上述三個入門例子中，Lax Milgram 定理告訴我們，他的解是存在唯一的。而對於一般的特徵值問題形式的方程，你再去驗證 Lax Milgram 的條件，你會發現，不一定滿足，這時候特徵值問題對應的解就不一定是存在唯一的，可能不存在，存在也不一定唯一。比如，對於$-\Delta u+bu=0$，只有 $b$ 是非負的時候， Lax Milgram 定理才成立。負的拉普拉斯運算元是一個半正定的運算元，加上 0 本質邊界條件，可以理解成了具有正定性。正定的運算元，有正的特徵值，所以當這裡的$b$取到拉普拉斯運算元的特徵值的時候，這個問題有非零解。這就是所謂的特徵值問題。可見，Lax Milgram定理只是充分條件，而不是必要條件。

你也可以用代數方程的情況去類比推想方程的情況。對於$Ax=bx$，$A$ 半正定。除零解外，若 $b$ 是特徵值，有一解多解的情況。若 $b$ 是負的不是特徵值，無解。考慮 $b=0$，如果 $0$ 是特徵值，有線性解 $bu+c$，限定邊界條件，這個線性解就定下來了。如果 $0$ 不是特徵值，無解。$0$ 特徵值只能是單特徵值，也就是說 $b$ 取零，本質邊界條件下，解都是唯一的。這段話比較拗口，可不看。

亥姆霍茲方程

亥姆霍茲方程 (英語: Helmholtz equation) 是一個描述電磁波的格圓偏微分方程, 以德國物理學家亥姆霍茲的名字命名。其基本形式如下:
$$\left({\nabla }^2+k^2\right)A=0$$
其中 ${\nabla }^2$ 是拉音拉斯運算元, $k$ 是波數, $A$ 是振幅。
因為它和波動方程的關係，亥姆霍茲方程在物理學中電磁輻射、地震學和聲學等相關研究領域裡有著廣泛應用。

程式與結果

問題一

求解問題： $$\frac{d^{2}u}{d{x}^2}+ku(x)=l(x)$$ 邊界條件：
$$\{\begin{array}{l}u(0)=0\\ u(1)=0\end{array}$$ 定義 $$l(x)=-\sin (\pi x)$$ 那麼解析解為
$$u(x)=\frac{1}{{\pi }^2-k}\sin (\pi x)$$ 當 $k=0.5{\pi }^2$時，解析解和數值解如圖

當 $k={\pi }^2$ 時，

當 $k=1.0001{\pi }^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

問題二

試著求解以下的亥姆霍茲方程，
$${\left(a{u}^{\prime }\right)}^{\prime }+k^{2}u=f\text{ in }(0,1)$$
邊界條件是，
$$\{\begin{array}{l}u(0)=0\\ u(1)=0\end{array}$$
這裡的
$$a(x):=1+\frac{1}{2}\sin \frac{x}{\epsilon }$$

$\epsilon$ 是週期的，$f$ 是光滑的。

測試你的程式碼對於
$$\epsilon =\frac{1}{2^n},n=0,1,2,3,4$$
通常地，為了得到精確解，你也需要改變你的網格尺寸：
$$h\sim \epsilon /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