有限元方法基础入门教程（一维弹性问题的有限元程序）

有限元方法（冯康首次发现时称为基于变分原理的差分方法），是一种用于求解微分方程组或积分方程组数值解的数值技术。这一解法基于完全消除微分方程，即将微分方程转化为代数方程组（稳定情形）；或将偏微分方程（组）改写为常微分方程（组）的逼近，这样可以用标准的数值技术（例如欧拉法，龙格－库塔法等）求解。

有限元方法的发展开始于五十年代中后期使用在机身框架和结构分析上，并于六十年代通过斯图加特大学的John Argyris和柏克莱加州大学的Ray W. Clough在土木工程中的应用工作中积累经验。1965年，冯康将其成果发表成论文。

FEM在理解起来并不困难，下面通过一个简单的例子，结合code，在说明这样一个问题。
1、首先考虑下面这样一个偏微分问题，物理建模过程略：

2、对于这个方程，任意找一个函数乘方程两边并进行积分，再利用分布积分公式，得到一个带积分的等式。进而，选取一组基函数，将uu和vv线性表出后代入，可将问题变为一个求解一个线性方程组问题。过程如下：

3、进一步，考虑到固有边界条件（u(0)=u0u(0)=u_0），最后要对刚度矩阵和载荷向量进行校正，如下：

有限元就是这么地简单，假如你认真看上面图片的话，应该能体会到。

下面给出这个问题的matlab代码：

主程序，调用了求解方程的函数：

%clear all; close all;L = 3; numel = 6; u_0 = 0; sigma_0 = 1; etype = 'linear';fem1d(L,numel,u_0,sigma_0,etype);% Exact and numerical solution will be plotted. Exact solution for these
% parameters can be found in exactsolution.m

fem1d.m如下，给出了一个求解该问题的函数，相信注释比较清楚了，多余的东西就不说了。

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Scientific Computing : Computational Methods for Nonlinear
% Eigenvalue Problems, LSEC, Chinese Academy of Sciences, Beijing, China.
%============================================================================
%
% Finite Element Code for One-Dimensional Elasticity
% --------------------------------------------------
% Author    : Song Lu, LSEC
% E-mail    : lusongcool@163.com
%
%============================================================================
%
% MODEL 1D Problem (Strong Form)
% ==============================
% One-dimensional Matlab finite element program to solve the following model
% problem (one-dimetttttttttttttttt6nsional elasticity):
%
% (EAu')' + b = 0  in \Omega = (0,L)
% where E(x) is the Young's modulus, A(x) is the cross-sectional area, and
% b(x) is the body force per unit length. The functions for E, A and b are
% set in m-files youngsmodulus.m, area.m and bodyforce.m, respectively.
%
% with boundary conditions
%
% u(0)   = u_0 (essential boundary condition)
% E(L)u'(L) = sigma_0 (natural boundary condition)
%
% FE Approximation
% ================
% Galerkin method with linear finite elements
%
% ----------------------------------------
% USAGE: fem1d(L,numel,u_0,sigma_0,etype)
% ----------------------------------------
% L       : length of the bar
% numel   : number of elements (nodes are equi-spaced)
% u_0     : essential BC at x = 0
% sigma_0 : traction at x = L
% etype   : element type (only option 'linear' is currently available)
% Default : L = 3; numel = 6; u_0 = 0; sigma_0 = 1; etype = 'linear'
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%function fem1d(L,numel,u_0,sigma_0,etype)clc;       % clear command window
close all; % close all figures% default valuesif nargin < 5, help(mfilename), %展示注释L = 3; numel = 6; u_0 = 0; sigma_0 = 1; etype = 'linear'; fprintf('Executing with default settings \n \n');
end if L < 0, error('L must be positive'); end
if numel < 1, error('Number of elements must be a natural number (1,2,...)'); endnumnod = numel + 1; % no_of_nodes = no_of_elements + 1 (linear elements)%
% check if element type (etype) is valid
%
switch lower(etype)%选择元的类型，目前只有线性的。case 'linear'nel = 2;disp('------------------------------------');disp('FEM analysis with linear elements');disp('------------------------------------');case 'quadratic'nel = 3;disp('---------------------------------------');disp('FEM analysis with quadratic elements');disp('---------------------------------------');error('Quadratic element type not coded yet');otherwisefprintf('Element type = %s not defined\n',etype);error('Aborting!');
end%
% Gauss quadrature (number of quadrature points)
%
nsp = 2; % Gauss-rule (nsp-point rule is used)%高斯点的个数% nodal discretizationcoord = linspace(0,L,numnod)'; % `numnod' equally-spaced nodes between 0 and L
%节点的坐标% nodal connectivity array for elements in the FE mesh (numel x 2 matrix)connect = [1:(numnod-1); 2:numnod]';%连接矩阵，表示哪些节点表示的基函数是连接的。
% essential boundariesebcdof = [1];   % first node at x = 0 is constrained
ebcval = [u_0]; % prescribed value
%在第一个点的地方为u0% initialize coefficient matrix [K] and external force vector {f}bigk = sparse(numnod,numnod); % stiffness matrix (initialized to zero) for [K]
fext = sparse(numnod,1);      % force vector (initialized to zero) for {f}%
% element loop BEGIN
%for e = 1:numelke = elemstiff(e,nel,nsp,coord,connect);fe = elemforce(e,nel,nsp,coord,connect);sctr = connect(e,:);
%
% assemble element stiffness matrix
%bigk(sctr,sctr) = bigk(sctr,sctr) + ke;
%
% assemble element force vector
%fext(sctr) = fext(sctr) + fe;
end%
% element loop END
%%
% incorporate the traction boundary term in the force vector (refer to the
% weak form)
%xend = coord(numnod);
fext(numnod) = fext(numnod) +  area(xend)*sigma_0;%
% essential boundary conditions
%% modify rows in K and f to incorporate the essential boundary condition
%根据本质边界条件进行校正
for i = 1:length(ebcdof)n = ebcdof(i);for j = 1:numnod if (isempty(find(ebcdof == j))) % dof j is not found in ebcdof vectorfext(j) = fext(j) - bigk(j,n)*ebcval(i);endendbigk(n,:) = 0.0;bigk(:,n) = 0.0;bigk(n,n) = 1.0;fext(n) = ebcval(i);
end%
% solve for nodal unknowns
%u_coeff = inv(bigk)*fext;%
% determine derivative of approximation in each element : for linear elements
% the strain (u') is a constant in each element and is given by the difference
% of the nodal displacements divided by the length of the element
%%求数值解的导数，因为用直线逼近，每一段都是一个常数
for e = 1:numelhe = coord(connect(e,2)) - coord(connect(e,1));uprime(2*e-1) = (u_coeff(connect(e,2)) - u_coeff(connect(e,1)))/he;uprime(2*e) =  uprime(2*e-1);xp(2*e-1) = coord(connect(e,1));  xp(2*e) = coord(connect(e,2));
end%
% exact solution : see exactsolution.m
%%
% number of sampling point to compute the exact solution for plotting purposes
nsamp = 101;
%
xsamp = linspace(0,L,nsamp);
for j = 1:length(xsamp)[uex,uexder] = exactsolution(xsamp(j));uexact(j) = uex;uprimeexact(j) = uexder;
end%
% plot results
%if (numel > 20) plot(coord,u_coeff,'-',xsamp,uexact,'--');
elseplot(coord,u_coeff,'-',xsamp,uexact,'--',...coord,zeros(numnod,1),'mo','MarkerEdgeColor','k',...'MarkerFaceColor',[0.5 1 0.6],'MarkerSize',10);
end
h = legend('FE solution','Exact solution','location','best');
set(h,'fontsize',14);
title('Comparison of FE and Exact Solutions');
xt = get(gca,'XTickLabel'); set(gca,'XTickLabel',xt,'fontsize',12);
yt = get(gca,'XTickLabel'); set(gca,'XTickLabel',yt,'fontsize',12);
xlabel('$x$','Interpreter','latex','fontsize',16)
ylabel('$u^h , u$','Interpreter','latex','fontsize',16);  figure;
if (numel > 20) plot(xp,uprime,'-',xsamp,uprimeexact,'--');
elseplot(xp,uprime,'-',xsamp,uprimeexact,'--',...coord,zeros(numnod,1),'mo','MarkerEdgeColor','k',...'MarkerFaceColor',[0.5 1 0.6],'MarkerSize',10);
end
h = legend('FE derivative','Exact derivative','location','best');
set(h,'fontsize',14);
title('Comparison of FE and Exact Derivative Solutions');
title('Comparison of FE and Exact Derivative Solutions');
xt = get(gca,'XTickLabel'); set(gca,'XTickLabel',xt,'fontsize',12);
yt = get(gca,'XTickLabel'); set(gca,'XTickLabel',yt,'fontsize',12);
xlabel('$x$','Interpreter','latex','fontsize',16)
ylabel('$u_{,x}^h , u_{,_x}$','Interpreter','latex','fontsize',16);return
end

一些子函数：

**elemstiff.m**function [ke] = elemstiff(e,nel,nsp,coord,connect)
%
% Purpose
% =======
% Element stiffness matrix for one-dimensional fem
%
% Input
% -----
% e       : element number
% nel     : nodes per element
% nsp     : no of Gauss points
% coord   : coordinates of the nodes
% connect : nodal connectivity
%
% Output
% ------
% Element stiffness matrix ke (nel x nel matrix)
%%
% initialize element stiffness matrix to zero
%
ke    = zeros(nel,nel); nodes = connect(e,:);
xe    = coord(nodes); % x-coordinate of nodes
Le    = xe(nel) - xe(1);%connect只有两列，这里为什么要用nel，考虑到隔一个也有交叉的基函数
detj  = Le/2;%
% loop over Gauss points in xi-direction
%xi = gauss_points(nsp);%-1到1上的高斯点
weight = gauss_weights(nsp);%高斯权重for i = 1:nsp
%
% shape function and derivative of N (B-matrix)
%N = [ 0.5*(1 - xi(i))  0.5*(1 + xi(i))];%基函数在高斯点上的值B = 1/Le*[-1 1];   % detj = L/2 求基函数的在高斯点的导数，即斜率xg = N*xe; %表示xe区间上高斯点的坐标
%
% assemble element stiffness matrix ke
%Ax = area(xg);%界面面积一般和位置也有关Ex = youngsmodulus(xg);%杨氏模量和位置有关ke = ke + weight(i)*Ex*Ax*B'*B*detj;
endreturn
end**elemforce.m**function [fe] = elemforce(e,nel,nsp,coord,connect)
%
% Purpose
% =======
% Element force vector for one-dimensional fem
%
% Input
% -----
% e       : element number
% nel     : nodes per element
% nsp     : no of Gauss points
% coord   : coordinates of the nodes
% connect : nodal connectivity
%
% Output
% ------
% Element force vector fe (nel x 1 vector)
%%
% initialize element force vector to zero
%
fe = zeros(nel,1); nodes = connect(e,:);
xe = coord(nodes); % x-coordinate of nodes
Le  = xe(nel) - xe(1);
detj = Le/2;%
% loop over Gauss points in xi-direction
%xi = gauss_points(nsp);
weight = gauss_weights(nsp);for i = 1:nsp
%
% shape function N
%N = [ 0.5*(1 - xi(i))  0.5*(1 + xi(i)) ];xg = N*xe;%
% assemble element force vector ke
%bx = bodyforce(xg);fe = fe + weight(i)*bx*N'*detj;
endreturn
end**guass_points.m**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;elseerror('GAUSS_WEIGHTS - Fatal error! Illegal value of n.')endreturn
end**gauss_weights.m**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;elseerror('GAUSS_WEIGHTS - Fatal error! Illegal value of n.')endreturn
end**area.m**
function Ax = area(xg)
%
% Cross-sectional area
%Ax = 1;return
end
**bodyforce.m**
function bx = bodyforce(xg)
%
% body force
%bx =  1;return
end
**youngsmodulus.m**
function Ex = youngsmodulus(xg)
%
% Young's modulus
%Ex = 1;return
end**exactsolution.m**function [uexact,uexactder] = exactsolution(xg)
%
% exact solution for L = 3, b = 1, E = 1, A = 1, u_0 = 0, sigma_0 = 1
%uexact = 4.*xg - xg^2/2.;
uexactder = 4. - xg;return
end

需要注意的是，杨氏模量等量需要根据实际物理问题给出，这里设置为常数。另外，高斯点和高斯权重是用于求积分的，程序中直接给出了数据，实际的求法可以做正交多项式，取次数最高项求零点即可，详见我的另外一篇博客。
高斯积分方法

在解偏微分方程的过程中，主要的难点是如何构造一个方程来逼近原本研究的方程，并且该过程还需要保持数值稳定性。目前有许多处理的方法，他们各有利弊。当区域改变时（就像一个边界可变的固体），当需要的精确度在整个区域上变化，或者当解缺少光滑性时，有限元方法是在复杂区域（像汽车、船体结构、输油管道）上解偏微分方程的一个很好的选择。例如，在正面碰撞仿真时，有可能在”重要”区域（例如汽车的前部）增加预先设定的精确度并在车辆的末尾减少精度（如此可以减少仿真所需消耗）;另一个例子是模拟地球的气候模式，预先设定陆地部分的精确度高于广阔海洋部分的精确度是非常重要的。

考虑到方便，变分推导过程就手写了；数学上的表述往往越严格越难理解，为了通俗易懂，有些地方故意说得不是那么严格，甚至懂的人会有看法；为了好理解，尽量避开了符号和专业术语。

有限元方法基础入门教程（一维弹性问题的有限元程序）相关推荐

Python基础入门教程：Day21-30/Web前端概述
Python基础入门教程:Web前端概述说明:本文使用的部分插图来自 Jon Duckett 先生的*HTML and CSS: Design and Build Websites*一书,这是一本非 ...
Android基础入门教程——8.3.18 Canvas API详解(Part 3)Matrix和drawBitmapMash
Android基础入门教程--8.3.18 Canvas API详解(Part 3)Matrix和drawBitmapMash 标签(空格分隔): Android基础入门教程本节引言: 在Canva ...
2021最新HarmonyOS鸿蒙系统应用开发之基础入门教程到实战—持续更新（第三节：鸿蒙的技术特征）
第三节:鸿蒙的技术特征每篇内容都有视频讲解,可直接点击观看+关注,持续更新中 2021最新HarmonyOS鸿蒙系统应用开发之基础入门教程到实战-持续更新(第二节:鸿蒙OS系统分布式操作) 硬件互助 ...
超全面Python基础入门教程【十天课程】博客笔记汇总表
目录 1.学习地址 2.下载地址[视频.源码.课件.软件] 3.详细博客笔记 Day 01 Day 02 Day 03 Day 04 Day 05 Day 06 Day 07 Day 08 Day 0 ...
Mongodb最基础入门教程
Mongodb最基础入门教程如果想了解一下redis的入门教程,可以去看一下我的上一篇博客 Mongodb的安装大家可以参考一下其他博主的博客,这里我就不做介绍了.不过值得注意的是,在Linux版本 ...
android设置webview缓存目录,Android基础入门教程——7.5.5 WebView缓存问题
Android基础入门教程--7.5.5 WebView缓存问题 Android基础入门教程本节引言:现在很多门户类信息网站,比如虎嗅,ifanr,钛媒体等等的APP,简单点说是信息阅读类的APP, ...
Python新手基础入门教程：如何像面试官展示你的coding能力
Python新手基础入门教程:如何像面试官展示你的coding能力作为程序员,我们通常在面试的时候面试官会出一些题来让你解决,作为优秀的程秀员,除了需要具备解决问题的思路以外,代码的质量也很关键.因 ...
视频教程-Python零基础入门教程-Python
Python零基础入门教程从2012年从事互联网至今有7年软件编程经验,曾任职国内北京互联网公司,中南林业大学授课Python 现任逻辑教育Python课程负责人,精通Python语言,精通人工智能 ...
SEO零基础入门教程（外链的发布和软文编写）
seo的作用是众所周知的,对网站进行seo优化,可以给网站带来大量的搜索引擎流量.但是想要做好网站优化也有难度,尤其是对于seo新手来说,因为缺乏理论和实战,所以seo新手需要多加练习.那么具体seo ...

有限元方法基础入门教程（一维弹性问题的有限元程序）

下面给出这个问题的matlab代码：

有限元方法基础入门教程（一维弹性问题的有限元程序）相关推荐

最新文章

热门文章