偏微分方程的有限元法求解
- 格式:doc
- 大小:823.51 KB
- 文档页数:22
偏微分方程数值求解方法偏微分方程数值求解方法是使用计算机算法来近似求解偏微分方程的过程。
偏微分方程是描述物理现象和自然现象的主要工具,但大多数偏微分方程不能通过解析方式求解,因此需要使用数值方法进行近似求解。
常用的偏微分方程数值求解方法包括有限差分法、有限元法、谱方法、边界元法和逆时空方法等。
1. 有限差分法有限差分法是一种最简单的数值求解方法,它将偏微分方程中的导数离散化为差分的形式,然后通过有限差分公式求解。
在有限差分法中,将求解区域离散化为网格,然后在每个节点上求解方程,通过节点之间的连通关系建立系数矩阵,最终利用线性代数方法求解线性方程组。
2. 有限元法有限元法是一种广泛运用的数值求解方法,它将求解区域离散化为有限个子域,然后在每个子域内近似求解方程。
有限元法是一种基于变分原理的方法,通过将偏微分方程转化为变分问题,然后在有限维的函数空间中建立逼近函数,最终利用变分方法求解方程。
3. 谱方法谱方法是一种基于傅里叶变换的数值求解方法,它将求解域上的函数表示为傅里叶级数的形式,然后通过求解系数来近似求解方程。
谱方法具有高精度、高效率的优点,但对于非周期边界和奇异性问题可能不适用。
4. 边界元法边界元法是一种基于积分方程的数值求解方法,它将偏微分方程转化为边界积分方程,然后在求解区域表面上求解方程。
边界元法不需要离散化求解区域,仅需在求解区域表面上采集节点,并通过节点之间的关系建立系数矩阵。
5. 逆时空方法逆时空方法是一种利用观测数据反演偏微分方程的数值求解方法,它通过最优化算法将观测数据反演为偏微分方程的参数。
逆时空方法对模型假设和观测数据的噪声较为敏感,但可以应用于各种偏微分方程的求解。
偏微分方程的数值求解方法偏微分方程是描述自然现象的重要工具,例如描述热传导、电磁波传播、流体运动等。
然而大多数情况下,这些方程很难通过解析方式求解,因此需要数值求解方法。
本文将介绍偏微分方程的数值求解方法及其应用。
一、有限差分法有限差分法是一种常见的偏微分方程数值求解方法。
它将原本连续的区域离散化,将偏微分方程转化为差分方程。
例如对于一维热传导方程:$$\frac{\partial u}{\partial t} = \alpha\frac{\partial^2 u}{\partial x^2} $$其中 $u(x, t)$ 是温度,$\alpha$ 是热扩散系数。
我们可以选择将空间分成 $N$ 个网格,时间分成 $M$ 个步骤。
则有:$$u_i^{m+1} = u_i^m + \frac{\alpha\Delta t}{\Deltax^2}(u_{i+1}^m - 2u_i^m + u_{i-1}^m)$$其中 $u_i^m$ 表示在位置 $i\Delta x$,时间 $m\Delta t$ 时的温度值。
这是一个显式求解方程,可以直接按照时间步骤迭代计算。
不过由于它的误差可能会增长,因此需要小心选择时间步长和空间步长,以保证误差不会过大。
二、有限元法有限元法是一种更加通用的偏微分方程数值求解方法。
它将连续区域离散化成一些小段,称为单元。
然后针对每个单元,将其上的偏微分方程转化为局部插值函数的方程求解。
例如对于一维波动方程:$$\frac{\partial^2 u}{\partial t^2} = c^2 \frac{\partial^2 u}{\partialx^2}$$我们可以选择将空间分成 $N$ 个网格,用有限元方法将每个网格分成若干个单元。
则对于每个单元 $i$,我们可以得到一个局部插值函数 $u^i(x, t)$ 来近似解该单元上的偏微分方程。
这里不再赘述该函数的形式。
另外,我们还需要满足界面上的连续性和斜率匹配条件,以保证整体解是连续的。
galerkin有限元法
galerkin有限元法
Galerkin有限元法,也称为Galerkin有限体积法(FV),是一种数值解决偏微分方程的有限元方法,用于快速求解各种椭圆型方程的数值求解。
它把椭圆型方程分解成多个有限元,然后对每个有限元计算其权重,将所有有限元的权重加起来就是椭圆型方程的数值解。
在使用Galerkin有限元法来解决椭圆型方程时,首先要确定有限元的形状与大小,这将影响有限元法求解时的准确程度。
一般来说,有限元的形状可以是矩形、三角形或其他任意多边形,但大小是由实际情况决定的,需要根据椭圆型方程质量结构以及实际求解精度来确定。
确定有限元的形状与大小之后,就可以为每个有限元应用Galerkin有限元法,主要步骤如下:
1. 对每个有限元确定一个适当的坐标系,以便计算其权重;
2. 将系数函数投影到有限元上,并且确定每个有限元的质点分布情况;
3. 确定每个有限元的权重,并将所有有限元的权重加起来就是椭圆型方程的数值解。
Galerkin有限元法的优点是可以快速求解出准确的解,而且可以灵活应用于解决多种椭圆型方程。
但是它也有一定的缺点,比如假设有限元的形状和大小得不到充分考虑,那么计算精度可能会降低;另外,在计算权重时,需要考虑每个有限元上的局部梯度,如果选取
的有限元尺度过小,必须计算大量的梯度,从而增加计算难度。
有限元分析如何求解偏微分方程有限元分析如何求解偏微分方程「篇一」1、有限元法是近似求解连续场问题的数值方法。
2、有限元法将连续的求解域(离散),得到有限个单元,单元与单元之间用(结点相连。
3、从选择未知量的角度看,有限元法可分为三类(位移法力法混合法)。
4、以(结点位移)为基本未知量的求解方法称为位移量。
5、以(结点力)为基本未知量的求解方法称为力法。
7、直梁在外力作用下,横截面上的内力有(剪力)和(弯矩)两个。
8、平面刚架结构在外力作用下,横截面上的内力有(剪力)、(弯矩)、(轴力)。
9、进行直梁有限元分析,结点位移有(转角)、(挠度)。
12、弹性力学问题的方程个数有(15)个,未知量个数有(15)个。
13、弹性力学平面问题方程个数有(8),未知数(8)个。
15、几何方程是研究(应变)和(位移)关系的方程。
16、物理方程描述(应力)和(应变)关系的方程。
17、平衡方程反映(应力)和(位移)关系的方程。
18、把进过物体内任意一点各个(截面)上的应力状况叫做(该点)的应力状态。
19、形函数在单元结点上的值,具有本点为(1),他点为零的性质,并在三角形单元的后一结点上,三个形函数之和为(1)。
20、形函数是(三角形)单元内部坐标的(线性位移)函数,它反映了单元的(位移)状态。
21、结点编号时,同一单元相邻结点的(编号)尽量小。
25、单元刚度矩阵描述了(结点力)和(结点位移)之间的关系。
矩形单元边界上位移是(线性)变化的。
1、从选择未知量的角度来看,有限元法可分为三类,下面那种方法不属于其中(C)。
A、力法B、位移法C、应变法D、混合法2、下面对有限元法特点的叙述中,哪种说法是错误的(D)。
A、可以模拟各种几何形状负责的结构,得出其近似值。
B、解题步骤可以系统化,标准化。
C、容易处理非均匀连续介质,可以求解非线性问题。
D、需要适用于整个结构的插值函数。
3、几何方程研究的是(A)之间关系的方程式。
A、应变和位移B、应力和体力C、应力和位移D、应力和应变 4.物理方研究的是(D)之间关系的方程式。
数值模拟偏微分方程的三种方法介绍(有限差分方法、有限元方法、有限体积方法)I.三者简介有限差分方法(Finite Difference Methods)是数值模拟偏微分方程最早采用的方法,至今仍被广泛使用O该方法包括区域剖分和差商代替导数两个步骤。
首先将求解区域划分为差分网格,用有限个网格节点代替连续的求解区域。
其次,利用Taylor级数展开等方法将偏微分方程中的导数项在网格节点上用函数值的差商代替进行离散,从而建立以网格节点上的值为未知量的代数方程组。
该方法是一种直接将微分问题变为代数问题的近似数值解法,数学概念直观,表达简单,是发展较早且十分成熟的数值方法。
差商代替导数后的格式称为有限差分格式,从格式的精度来考虑,有一阶格式、二阶格式和高阶格式。
从差分的空间离散形式来考虑,有中心格式和迎风格式。
对于瞬态方程,考虑时间方向的离散,有显格式、隐格式、交替显隐格式等。
目前常见的差分格式,主要是以上几种格式的组合,不同的组合构成不同的差分格式。
差分方法主要适用于结构网格,网格的大小一般根据问题模型和Courant稳定条件来决定。
有限元方法(Finite Element Methods)的基础是虚位移原理和分片多项式插值。
该方法的构造过程包括以下三个步骤。
首先,利用虚位移原理得到偏微分方程的弱形式,将计算区域划分为有限个互不重叠的单元(三角形、四边形、四面体、六面体等),在每个单元上选择合适的节点作为求解函数的插值点,将偏微分方程中的变量改写成由各变量或其导数的节点值与所选用的分片插值基函数组成的线性表达式,得到微分方程的离散形式。
利用插值函数的局部支集性质及数值积分可以得到未知量的代数方程组。
有限元方法有较完善的理论基础,具有求解区域灵活(复杂区域)、单元类型灵活(适于结构网格和非结构网格)、程序代码通用(数值模拟软件多数基于有限元方法)等特点。
有限元方法最早应用于结构力学,随着计算机的发展已经渗透到计算物理、流体力学与电磁学等各个数值模拟领域。
数值计算中的偏微分方程解法偏微分方程在科学、工程和金融等领域都有广泛的应用。
在现实生活中,许多问题都涉及到偏微分方程的解法,比如天气预报、机器学习和金融衍生品定价等。
然而,解析解并不总是可行的,因此需要数值计算方法来解决这些问题。
在本文中,我们将探讨数值计算中的偏微分方程解法。
一、有限差分法有限差分法是偏微分方程数值解法中最基本的方法之一。
该方法通过将偏微分方程中的导数用差分近似公式表示出来,然后建立一个离散的空间和时间网格。
在网格上求解方程,得到数值解。
例如,考虑一个二维热传导方程:$$ \frac{\partial u}{\partial t}= \alpha \left( \frac{\partial ^2u}{\partial x^2} +\frac{\partial ^2 u}{\partial y^2} \right) $$其中,$u(x,y,t)$是温度分布,$\alpha$是热传导系数。
我们可以将该方程在空间上进行离散化,用差分近似公式表示出导数。
以二阶中心差分为例,有:$$ \frac{\partial ^2 u}{\partial x^2} \approx \frac{u_{i+1,j}-2u_{i,j}+u_{i-1,j}}{\Delta x^2} $$$$ \frac{\partial ^2 u}{\partial y^2} \approx \frac{u_{i,j+1}-2u_{i,j}+u_{i,j-1}}{\Delta y^2} $$其中,$u_{i,j}$表示网格点$(i,j)$处的温度。
同样地,时间上也进行离散化,用前向差分公式表示导数,即:$$ \frac{\partial u}{\partial t} \approx \frac{u_{i,j}^{n+1}-u_{i,j}^n}{\Delta t} $$将上述离散化的结果代入方程中,可以得到:$$ \frac{u_{i,j}^{n+1}-u_{i,j}^n}{\Delta t}= \alpha\left( \frac{u_{i+1,j}^n-2u_{i,j}^n+u_{i-1,j}^n}{\Delta x^2}+\frac{u_{i,j+1}^n-2u_{i,j}^n+u_{i,j-1}^n}{\Delta y^2} \right) $$整理得到:$$ u_{i,j}^{n+1}= u_{i,j}^n+ \frac{\alpha \Delta t}{\Delta x^2} (u_{i+1,j}^n-2u_{i,j}^n+u_{i-1,j}^n)+ \frac{\alpha \Delta t}{\Delta y^2} (u_{i,j+1}^n-2u_{i,j}^n+u_{i,j-1}^n) $$这样,我们就可以用迭代法求解上述方程,得到网格上的温度分布。
16.901讲义笔记一维有限%首先,我们考虑•个比上一节稍微复杂点的问题; 豎二f(X),卫冲,V(O) = O.V(L)=O在这里,f(X)是)C的般函数,我们来看•个特别的情形:f(x)=x(L-x),此时,方程的梏确解如F:有限元方法利用加权残差的方法■其中:(1)设va)=£«Ma), v()()是我们对v(x)的近似,省为未知常数9 V|(x)是用户选择的歯数,即形状朗数:(2)定义N个加权残差LRj = p^(x)R(V)dx • j = l-> N to其中,RV)二器・f为绒差凹⑴足“用户”选择的加权函数,即权函数:(3)令加权残并为冬•町以确定⑷的值,即求耳使得对所fi 1=I->N, Rj=Oe令限元方法( )是加权残若法的一种,下血看看我们是如何用它来解决问题的。
一维有限元方法有限元方法(〉扌野个连续区域离散化-系列小单尤,这些单元与有限差分法()或有限体积法()产牛的网格完全相同,而佼之前两者主耍的优点在于:能够容易地把握单元的变化范囤。
对于我们讨论的一维问题,可以将区域(数轴〉离散化为如下图所示:这里,叫三单•元的个数。
我们还会用別下血i些定义:个三角划分;尽管令限元法对于一维,二维,三维甚至高细问题都是仃效的,们我们还是要谈及区域离散化的一种方浓,即三角划分。
4 T定义为第I个单元所在的区域。
对于_维问题,这表明,TS-个满足片心的X的集合。
接卜来耍确定的是毎个单兀该用什么样的函数,典型的函数形式就是用从一个单元到卜一个单兀保持解连续的多项式。
例如:一个线性有限元如卜團;i示:在毎个单元内的函数是线形的,在毎两个单元的交点处足连续的。
对于专门诜择的满足线件变化的形状函数,右估计残差时有一个很明显的问题:回忆前曲的内容,RV)二器一f,它在一个单冗里等于什么呢?因为函数是线性的,所以器=0,则有:R(V)=f ,即R(V)与无关。
冋时,满足线性变化的形状函数似乎也是一个好的近似,我们举-个例子来说明。
有限元法与偏微分方程的数值解法在现代科学技术中,物理和工程问题通常涉及到方程的解析解。
然而,有很多复杂的问题,没有精确的解析解。
在这些情况下,我们可以使用数值方法来解决问题。
其中,有限元法(Finite Element Method,FEM)被广泛应用于求解偏微分方程(Partial Differential Equation,PDE)的数值解法。
有限元法是一种数值解法,用于解决连续介质(如固体、液体和气体)的差分方程。
它通常涉及将整个计算域分成许多小区域,称为有限元。
这些有限元被视为形状简单的几何单元(如三角形、四边形、六边形等),并且为每个元素分配了未知值。
在有限元方法中,偏微分方程被转换为一个离散方程,其中未知数在局部有限元中定义。
该方法通常涉及将初始有限元网格粗略地分配到整个计算区域,以构建数值解的近似值。
我们可以使用数学方法,如高斯消元法或迭代方法,来求解这个离散的线性系统。
有限元方法在许多领域中发挥着重要作用,包括结构力学、流体力学、电磁学、信息学和生物工程等。
它可以用于求解几乎所有类型的PDE,例如:椭圆、双曲和抛物型等。
在有限元方法中,解取决于网格的精度。
对于较小的网格,精度较高,但计算时间较长;反之亦然。
因此,在选择网格时需要进行权衡。
此外,一个好的网格应该是稳定的,能够保证数值解的收敛性和精度。
一些常见的有限元方法包括:显式和隐式欧拉方法、二阶Runge-Kutta 方法和高阶方法等。
这些方法主要涉及将初始条件和边界条件应用到整个计算区域。
作为一种广泛使用的数值解法,有限元法已经成为许多计算机辅助工程计算软件的主要工具,例如有限元分析软件 ANSYS 等。
此外,计算机的性能提高了许多,使得我们能够处理更多的网格和更大的计算域。
结论有限元法是一种强大的数值解法,可用于求解广泛的物理和工程问题。
然而,对于不同的应用,有不同的适用条件和精度要求。
因此,在设计计算方案之前,需要进行仔细的分析和权衡,以确保最终的数值解具有良好的收敛性和精度。
二阶偏微分方程的 Matlab有限元法求解摘要:本文基于偏微分方程有限元法求解原理,运用Matlab中的偏微分方程工具箱(PDE Toolbox)对三类典型的二阶偏微分方程:椭圆型方程、双曲线型方程和抛物线型方程算例进行求解,为求解偏微分方程的提供参考。
关键词:偏微分方程,有限元,Matlab偏微分方程工具箱0引言偏微分方程定解问题是描述许多自然现象或工程问题的最重要的数学模型,应用非常广泛[1]。
解析法只能求解非常简单的偏微分方程,远远不能满足科学研究和工程实际的需要。
随着计算机技术和科学计算的迅速发展,数值解法成为求解偏微分方程的重要工具[2-3]。
数值解法将连续问题离散化,最后将偏微分方程化成线性代数方程组。
根据离散化方法不同,偏微分方程数值解法主要有差分法和有限元法。
有限元法是分片定义试函数与变分原理相结合的产物。
它能适应各种形状的区域,且通用性强,现已成为求解偏微分方程定解问题的一种有效数值方法[4]。
本文首先简述了偏微分方程有限元法原理,然后,对Matlab中的偏微分方程工具箱(Partial Differential Equations Toolbox)的功能和求解思路进行了阐述[5-6],最后,给出了用PDE Toolbox求解椭圆方程、、双曲线方程和抛物线方程的计算实例。
1偏微分方程有限元法原理偏微分方程有限元法的基本思想是将实际上连续的整个求解域进行离散化处理,即用一些假想的面或线将求解域分割为一系列的单元,各个单元之间仅在有限个节点处相互连接。
取未知函数的节点值作为基本未知量,在每个单元上选取一个近似的插值函数表示单元中场函数的分布规律。
利用变分原理来获得单元的刚度方程,然后按一定的规则把所有单元的刚度方程组集合起来,经适当的边界条件处理,便得到整个系统的总体方程组。
这样,偏微分方程便转化为一组常微分方程。
最后,求解总体方程组,得到节点值和用插值函数确定整个求解域上的场函数。
偏微分方程的数值解法偏微分方程(Partial Differential Equation, PDE)是数学和物理学中的重要概念,广泛应用于工程、科学和其他领域。
在很多情况下,准确解析解并不容易获得,因此需要利用数值方法求解偏微分方程。
本文将介绍几种常用的数值解法。
1. 有限差分法(Finite Difference Method)有限差分法是最常见和经典的数值解法之一。
基本思想是将偏微分方程在求解域上进行离散化,然后用差分近似代替微分运算。
通过求解差分方程组得到数值解。
有限差分法适用于边界条件简单且求解域规则的问题。
2. 有限元法(Finite Element Method)有限元法是适用于不规则边界条件和求解域的数值解法。
将求解域划分为多个小区域,并在每个小区域内选择适当的形状函数。
通过将整个域看作这些小区域的组合来逼近原始方程,从而得到一个线性代数方程组。
有限元法具有较高的灵活性和适用性。
3. 有限体积法(Finite Volume Method)有限体积法是一种较新的数值解法,特别适用于物理量守恒问题。
它通过将求解域划分为多个控制体积,并在每个体积内计算守恒量的通量,来建立离散的方程。
通过求解这个方程组得到数值解。
有限体积法在处理守恒律方程和非结构化网格上有很大优势。
4. 局部网格法(Local Grid Method)局部网格法是一种多尺度分析方法,适用于具有高频振荡解的偏微分方程。
它将计算域划分为全局细网格和局部粗网格。
在全局细网格上进行计算,并在局部粗网格上进行局部评估。
通过对不同尺度的解进行耦合,得到更精确的数值解。
5. 谱方法(Spectral Method)谱方法是一种基于傅里叶级数展开的高精度数值解法。
通过选择适当的基函数来近似求解函数,将偏微分方程转化为代数方程。
谱方法在处理平滑解和周期性边界条件的问题上表现出色,但对于非平滑解和不连续解的情况可能会遇到困难。
6. 迭代法(Iterative Method)迭代法是一种通过多次迭代来逐步逼近精确解的求解方法。
在科学技术各领域中,有很多问题都可以归结为偏微分方程问题。
在物理专业的力学、热学、电学、光学、近代物理课程中都可遇见偏微分方程。
偏微分方程,再加上边界条件、初始条件构成的数学模型,只有在很特殊情况下才可求得解析解。
随着计算机技术的发展,采用数值计算方法,可以得到其数值解。
偏微分方程基本形式而以上的偏微分方程都能利用PDE工具箱求解。
PDE工具箱PDE工具箱的使用步骤体现了有限元法求解问题的基本思路,包括如下基本步骤:1) 建立几何模型2) 定义边界条件3) 定义PDE类型和PDE系数4) 三角形网格划分5) 有限元求解6) 解的图形表达以上步骤充分体现在PDE工具箱的菜单栏和工具栏顺序上,如下具体实现如下。
打开工具箱输入pdetool可以打开偏微分方程求解工具箱,如下首先需要选择应用模式,工具箱根据实际问题的不同提供了很多应用模式,用户可以基于适当的模式进行建模和分析。
在Options菜单的Application菜单项下可以做选择,如下或者直接在工具栏上选择,如下列表框中各应用模式的意义为:① Generic Scalar:一般标量模式(为默认选项)。
② Generic System:一般系统模式。
③ Structural Mech.,Plane Stress:结构力学平面应力。
④ Structural Mech.,Plane Strain:结构力学平面应变。
⑤ Electrostatics:静电学。
⑥ Magnetostatics:电磁学。
⑦ Ac Power Electromagnetics:交流电电磁学。
⑧ Conductive Media DC:直流导电介质。
⑨ Heat Tranfer:热传导。
⑩ Diffusion:扩散。
可以根据自己的具体问题做相应的选择,这里要求解偏微分方程,故使用默认值。
此外,对于其他具体的工程应用模式,此工具箱已经发展到了Comsol Multiphysics软件,它提供了更强大的建模、求解功能。
%对于d2u/dx2=f的FEM解算器,其中f=x*(1-x)%%边界条件u(0)=0, u(1)=0.%精确解用以比对xx=linspace(0,1,101);%产生0-1之间的均分指令,101为元素个数uex=(1/6).*xx.^3-(1/12).*xx.^4-(1/12).*xx;%对力项设置高斯点的数目NGf=2;if (NGf==2)xiGf=[-1/sqrt(3);1/sqrt(3)];%ξ1、ξ2的值aGf=[1 1];else,NGf=1;xiGf=[0.0];aGf=[2.0];end%单元数目Ne=5;%建立网格节点x=linspace(0,1,Ne+1);%零刚性矩阵K=zeros(Ne+1,Ne+1);b=zeros(Ne+1,1);%对所有单元循环计算刚性和残差for ii=1:Ne,kn1=ii;kn2=ii+1;x1=x(kn1);x2=x(kn2);dx=x2-x1;%每一个单元的长度dxidx=2/dx;%dξ/dxdxdxi=1/dxidx;%dx/dξdN1dxi=-1/2;%dζ1/dξdN2dxi=1/2;%dζ2/dξdN1dx=dN1dxi*dxidx;%-1/(x j-x j-1)dN2dx=dN2dxi*dxidx;%1/(x j-x j-1)K(kn1,kn1)=K(kn1,kn1)-2*dN1dx*dN1dx*dxdxi;%Rj的第二项K(kn1,kn2)=K(kn1,kn2)-2*dN1dx*dN2dx*dxdxi;K(kn2,kn1)=K(kn2,kn1)-2*dN2dx*dN1dx*dxdxi;K(kn2,kn2)=K(kn2,kn2)-2*dN2dx*dN2dx*dxdxi;%用高斯积分估计力项的积分for nn=1:NGf%NGf=2xiG=xiGf(nn);%得到高斯点的ξN1=0.5*(1-xiG);%求N1和N2(即在xiG的权重/插值) 形状函数在ξ的值N2=0.5*(1+xiG);%ζ值fG=xiG*(1-xiG);%对ξ点求fgG1=N1*fG*dxdxi;%在节点处估计权函数在高斯点的被积函数gG2=N2*fG*dxdxi;%估计是个积分值b(kn1)=b(kn1)+aGf(nn)*gG1;% aGf为1b(kn2)=b(kn1)+aGf(nn)*gG2;endend%在x=0处设置Dirichlet条件kn1=1;K(kn1,:)=zeros(size(1,Ne+1));K(kn1,kn1)=1;b(kn1)=0;%在x=1处设置Dirichlet条件kn1=1;K(kn1,:)=zeros(size(1,Ne+1));K(kn1,kn1)=1;b(kn1)=0;%求解方程v=K\b;%v为Kx=b的解plot(x,v,'*-');%画图并比较hold on;plot(xx,uex);hold off;xlabel('x');ylabel('u');。
%对于d2u/dx2=f的FEM解算器,其中f=x*(1-x)
%
%边界条件u(0)=0, u(1)=0.
%精确解用以比对
xx=linspace(0,1,101);%产生0-1之间的均分指令,101为元素个数
uex=(1/6).*xx.^3-(1/12).*xx.^4-(1/12).*xx;
%对力项设置高斯点的数目
NGf=2;
if (NGf==2)
xiGf=[-1/sqrt(3);1/sqrt(3)];%ξ1、ξ2的值
aGf=[1 1];
else,
NGf=1;
xiGf=[0.0];
aGf=[2.0];
end
%单元数目
Ne=5;
%建立网格节点
x=linspace(0,1,Ne+1);
%零刚性矩阵
K=zeros(Ne+1,Ne+1);
b=zeros(Ne+1,1);
%对所有单元循环计算刚性和残差
for ii=1:Ne,
kn1=ii;
kn2=ii+1;
x1=x(kn1);
x2=x(kn2);
dx=x2-x1;%每一个单元的长度
dxidx=2/dx;%dξ/dx
dxdxi=1/dxidx;%dx/dξ
dN1dxi=-1/2;%dζ1/dξ
dN2dxi=1/2;%dζ2/dξ
dN1dx=dN1dxi*dxidx;%-1/(x j-x j-1)
dN2dx=dN2dxi*dxidx;%1/(x j-x j-1)
K(kn1,kn1)=K(kn1,kn1)-2*dN1dx*dN1dx*dxdxi;%Rj的第二项
K(kn1,kn2)=K(kn1,kn2)-2*dN1dx*dN2dx*dxdxi;
K(kn2,kn1)=K(kn2,kn1)-2*dN2dx*dN1dx*dxdxi;
K(kn2,kn2)=K(kn2,kn2)-2*dN2dx*dN2dx*dxdxi;
%用高斯积分估计力项的积分
for nn=1:NGf%NGf=2
xiG=xiGf(nn);%得到高斯点的ξ
N1=0.5*(1-xiG);%求N1和N2(即在xiG的权重/插值) 形状函数在ξ的值
N2=0.5*(1+xiG);%ζ值
fG=xiG*(1-xiG);%对ξ点求f
gG1=N1*fG*dxdxi;%在节点处估计权函数在高斯点的被积函数gG2=N2*fG*dxdxi;%估计是个积分值
b(kn1)=b(kn1)+aGf(nn)*gG1;% aGf为1
b(kn2)=b(kn1)+aGf(nn)*gG2;
end
end
%在x=0处设置Dirichlet条件
kn1=1;
K(kn1,:)=zeros(size(1,Ne+1));
K(kn1,kn1)=1;
b(kn1)=0;
%在x=1处设置Dirichlet条件
kn1=1;
K(kn1,:)=zeros(size(1,Ne+1));
K(kn1,kn1)=1;
b(kn1)=0;
%求解方程
v=K\b;%v为Kx=b的解
plot(x,v,'*-');%画图并比较
hold on;
plot(xx,uex);
hold off;
xlabel('x');
ylabel('u');。