有限单元法学习报告
在对力学问题分析求解过程中，方法可以概括为两种方法，一种为解析法，对具体问题具体分析，通过一定的推导用具体的表达式获得解答，由于实际工程中结构物的复杂性，此方法在处理工程问题是十分困难的；另一种是数值法，有限元法是其中一种方法，其数学逻辑严谨，物理概念清晰，又采用矩阵形式表达基本公式，便于计算机编程，因此在工程问题中获得广泛的应用。
有限元法基本原理是，将复杂的连续体划分为简单的单元体；将无限自由度问题化为有限自由度问题，因为单元体个数是有限的；将偏微分方程求解问题化为有限个代数方程组的求解问题。通常以位移为基本未知量，通过虚功原理和最小势能原理来求解。
基本思想是先化整为零，即离散化整体结构，把整体结构看作是由若干个通过结点相连的单元体组成的整体；再积零为整，通过结点的平衡来建立代数方程组，最后计算出结果。
我将采用最简单的三结点三角形为基本单元体，解决弹性力学中的平面问题为例，解释有限单元法的基本原理、演示数值计算过程和一般性应用结论。
一、离散化
解决平面问题时，主要单元类型包括三角形单元（三结点、六结点）和四边形单元（四结点矩形、四结点四边形、八结
点四边形）等。选用不同的单元会有不同的精度，划分的单元数越多，精度越高，但计算量也会越大。因此在边界曲折，应力集中处单元的尺寸要小些，但最大与最小单元的尺寸倍数不宜过大。在集中力作用点及分布力突变的点宜选为结点，不同厚度，不同材料不能划分在同一单元中。三角形单元以内角接近60°为最好。充分利用对称性与反对称性。二、单元分析
将一个单元上的所有未知量用结点位移表示，并将分布在单元上的外力等效到结点上。
1、位移函数选取：
根据有限元法的基本思路，将连续体离散为有限的单元集合后，此时单元体满足连续性、均匀性、各向同性、完全线弹性假设。单元与单元之间通过结点连接并传递力，位移法（应用最广）以结点位移δi=(u i v i)T为基本未知量，以离散位移场代替连续位移场。单元体内的位移变化可以用位移函数（位移模式）来表示，因为有限元分析所得结果是近似结果，为了保证计算精度和收敛性，x位移函数应尽可能反应物体中的真实位移，即满足完备性和连续性的要求：
①位移模式必须能反映单元的刚体位移。
②位移模式必须能反映单元的常量应变。
③位移模式应尽可能反应位移的连续性。
设三角形单元三个结点编号为i、j、m。平面三角形单元位
移函数选取为
u=α1+α2x+α3y v=α4+α5x+α6y
可以写成00u u y
v v y ωω=-??=+?
的形式，00u v 、反映了单元的刚体平动，ω反
映了单元的刚体转动，满足完备性和连续性的要求①。 采用插值法由单元结点位移列阵δe
=()
i
i j j m
m u v u v u v T
计
算α1、α2、α3、α4、α5、α6.，求出位移d=[u （x ，y), v （x ，y ）]。6个未知量,6个代数方程,得d e
=N δe
d
e
=u v ?? ???=0000
0i j m i
j
m N N N N N N ?? ???()
i
i j j m
m u v u v u v T
式中N i =(a i +b i x+c i y)/2A ，a i =j j m
m
x y x y b i = -
11i
m
y y c i =
11j
m
x x （i 、j 、m 轮换）A 为三角形面积，为避免A<0，i 、j 、m 按逆时针排列。N 为形函数矩阵，形函数Ni 的性质有：
①N i （x i ，y i ）=1 N i （x j ，y j ）=0 Ni （xm ，ym ）=0 ②N i （x ，y ）+N j （x ，y ）+N m （x ，y ）=1可推出三个形函数中，两个是独立的，反映了刚体平移。
令z=Ni ，在直接坐标系中画出Ni 、Nj 、Nm 的函数图形是以Ni （xi ，yi ）=1为高的四面体，所以结点位移影响单元的位移场，单元的位移场是线性分布的，相邻单元在公共边上的位移是连续的，单元相邻边的位移只取决于单元相邻公共边上的结点而与其他结点无关，无论以哪个单元计算相邻边的位移，结果一定相同。
形函数N i e
决定了单元内的位移模式，反映了i 结点位移对单元内任意点位移的贡献率。
2、根据几何方程用单元结点位移表示单元应变：
()
00010002i
j m i
j m i i j j m
m i i
j
j
m
m u x B b b v c c c u v u v u v y A c
b c b c b u v y x ε??
? ?
? ??? ?? ?== ? ?? ? ???
?
??+ ????
?T
e e B εδ=
B 为几何矩阵
B 可写为分块矩阵
B =（B i B j B m ）T
,B i =00
i
i i
i b c c b ??
?
? ???
，B 内所有元素
与x ，y 无关，所以该单元内应变是常量，反映单元的常量应变，满足完备性和连续性的要求 ，这是一种常应变单元。 3、根据物理方程用单元结点位移表示单元应力：
e e D σε= 210101100
2E D μ
μμμ?
? ?
?
= ?- ?
- ??
?
D 为弹性矩阵
e e e e D DB S σεδδ=== S 为应力矩阵
S=DB 中，每一个元素都是常数，所以e σ的每一个分量与单元内x ，y 位置无关，这是一种常应力单元。
因为在三结点三角形单元中，位移函数中含有坐标的一次项，其误差为()2o x ?，而应力、
应变是常量，其误差为()o x ?，比位移精度低。 4、根据虚功原理用单元结点位移表示单元结点力
单元在结点处受力，单元会发生变形，因此单元在结点处所受到的力与单元结点位移肯定有关系。单元间通过结点的相互作用成为整体，因此每一单元的受力——位移关系找出来，整体的受力——位移关系也就出来了。 记单元节点力为()
e i
j
m F F F F =T
，单元结点虚位移为
()
*e
***i j m =T
δδδδ
单元内应力为()
e x
y xy σσστ=T
， 单元内虚应变
()
****x y
xy εεεγ=T
根据虚功原理，()
()**T
T
e e e A
F dxdy t δεσ=???，可得
e T e A
F B DBdxdy t δ=????
因为B 、D 中元素都是常数，e T e e F B DBtA K δδ==，K=B T
DB tA 为单元刚度矩阵。 K 为6行6
列矩阵可写为()T i T j i j
m T m B K tA B D B B B B ?? ?== ? ???
ii
ij im ji
j j jm mi
mj
mm k k k k k k k k k ?? ? ? ???
，xx
xy ij ij ij i j yx
yy ij
ij k k k B DB At At k k ??== ? ???
，xy
ij k 表示j 结点处发生y 方向的单位位
移时所引起的i 结点处x 方向的结点力。不同类型不同形式的单元，只有弹性矩阵D 和几何矩阵B 不同，计算子块矩阵的公式相同，平面问题中，影响刚度矩阵K 的只有几何矩阵B 。
K 的性质有：
K 中每个元素表示个单元结点沿坐标方向发生单位位移时所引起的结点力。
②K 为对称矩阵。
③单元做刚体位移时，单元内不产生应变应力，结点力为0，所以K 中每行每列元素之和为0，所以0K
=，所以只根据
e e F K δ=无法求得唯一解。
5、根据虚功等效原则计算等效结点力
根据有限元的基本方法，单元内任意点的位移、应变、应力等最终都要用结点位移来表示，所以作用在物体上的外力也要用结点位移表示。为了计算等效结点力，在任意的虚位移上，使原载荷与等效载荷虚功相等。
设外力为p f ，结点虚位移为*e δ，则任意点虚位移为**e e d N δ=，等效节点载荷为e L F ,有
*eT *eT e p L d f t F δ= e T L p F N f t =（集中力）
同理得e T L S F N f ds t =??（面力），e L A F Nfdxdy t =???（体力）。 三、整体分析
将结构的所有单元通过结点连接起来，形成一个整体的离散结构以代替实际的连续体，以形成以结点位移为未知量的整体结构的有限元代数方程组，最后求得结点位移。 对结点受力分析：结点受到与之相关的单元给它的反作用力和外载荷的等效结点力，这两组力坐标轴方向相反，所以应该相等，即i Li e
e
F F =∑∑，设有n 个结点，每个结点建立两个方
向的方程，不考虑外界约束时，共2n 个方程，2n 个未知量（,,...ix iy jx δδδ），为了建立这个代数方程组，建立整个弹性体的
结点力和结点位移的关系式L K F δ=，K （2n ×2n ）为整理刚度矩阵，δ为整体结点位移列阵，F L 整体结点载荷列阵。 为了求整体刚度矩阵，要找到它与已求得的单元刚度矩阵的关系，在整体中对结点编码，设整体刚度矩阵中某元素为Kij ，意为j 个结点在x 或y 方向发生位移引起i 个结点x 或y 方向的结点力，找到同时用到i 与j 结点的单元，并用与之对应的单元刚度矩阵中的元素ksm 相加得到Kij ，整体刚度矩阵也是奇异矩阵，必须考虑边界约束条件，消除K 的奇异性，才能求解结点位移。再由单位结点等效载荷得到整体结点载荷列阵F L 。这样K 、F L 已知，求解代数方程，解出整体结点位移列阵δ，得到相应的单元结点位移δe
。 δe
得到了，相应的d e
、σe
、εe
等就得到了。
【 数值计算方法试题一 一、 填空题（每空1分，共17分） 1、如果用二分法求方程043=-+x x 在区间]2,1[内的根精确到三位小数，需对分（ ）次。 2、迭代格式)2(2 1-+=+k k k x x x α局部收敛的充分条件是α取值在（ ）。 3、已知?????≤≤+-+-+-≤≤=31)1()1()1(211 0)(2 33x c x b x a x x x x S 是三次样条函数， 则 a =( )， b =（ ）， c =（ ）。 4、)(,),(),(10x l x l x l n 是以整数点n x x x ,,,10 为节点的Lagrange 插值基函数，则 ∑== n k k x l 0)(( )， ∑== n k k j k x l x 0 )(( )，当2≥n 时 = ++∑=)()3(20 4x l x x k k n k k ( )。 ； 5、设1326)(2 47+++=x x x x f 和节点,,2,1,0,2/ ==k k x k 则=],,,[10n x x x f 和=?07 f 。 6、5个节点的牛顿-柯特斯求积公式的代数精度为 ，5个节点的求积公式最高代数精度为 。 7、{}∞ =0)(k k x ?是区间]1,0[上权函数x x =)(ρ的最高项系数为1的正交多项式族，其中1)(0=x ?，则?= 1 4)(dx x x ? 。 8、给定方程组?? ?=+-=-2211 21b x ax b ax x ，a 为实数，当a 满足 ，且20<<ω时，SOR 迭代法收敛。 9、解初值问题 00 (,)()y f x y y x y '=?? =?的改进欧拉法 ??? ??++=+=++++)],(),([2),(] 0[111] 0[1n n n n n n n n n n y x f y x f h y y y x hf y y 是 阶方法。
有限单元法学习报告 在对力学问题分析求解过程中，方法可以概括为两种方法，一种为解析法，对具体问题具体分析，通过一定的推导用具体的表达式获得解答，由于实际工程中结构物的复杂性，此方法在处理工程问题是十分困难的；另一种是数值法，有限元法是其中一种方法，其数学逻辑严谨，物理概念清晰，又采用矩阵形式表达基本公式，便于计算机编程，因此在工程问题中获得广泛的应用。 有限元法基本原理是，将复杂的连续体划分为简单的单元体；将无限自由度问题化为有限自由度问题，因为单元体个数是有限的；将偏微分方程求解问题化为有限个代数方程组的求解问题。通常以位移为基本未知量，通过虚功原理和最小势能原理来求解。 基本思想是先化整为零，即离散化整体结构，把整体结构看作是由若干个通过结点相连的单元体组成的整体；再积零为整，通过结点的平衡来建立代数方程组，最后计算出结果。我将采用最简单的三结点三角形为基本单元体，解决弹性力学中的平面问题为例，解释有限单元法的基本原理、演示数值计算过程和一般性应用结论。 一、离散化 解决平面问题时，主要单元类型包括三角形单元（三结点、六结点）和四边形单元（四结点矩形、四结点四边形、八结点四边形）等。选用不同的单元会有不同的精度，划分的单元数越多，精度越高，但计算量也会越大。因此在边界曲折，应力集中处单元的尺寸要小些，但最大与最小单元的尺寸倍数不宜过大。在集中力作用点及分布力突变的点宜选为结点，不同厚度，不同材料不能划分在同一单元中。三角形单元以内角接近60°为最好。充分利用对称性与反对称性。 二、单元分析 将一个单元上的所有未知量用结点位移表示，并将分布在单元上的外力等效到结点上。 1、位移函数选取： 根据有限元法的基本思路，将连续体离散为有限的单元集合后，此时单元体满足连续性、均匀性、各向同性、完全线弹性假设。单元与单元之间通过结点连接并传递力，位移法（应用最广）以结点位移δi=(u i v i)T为基本未知量，以离散位移场代替连续位移场。单元体内的位移变化可以用位移函数（位移模式）来表示，因为有限元分析所得结果是近似结果，为了保证计算精度和收敛性，x位移函数应尽可能反应物体中的真实位移，即满足完备性和连续性的要求：
习 题 四 解 答 1、设010,1x x ==，写出()x f x e -=的一次插值多项式1()L x ，并估计插值误差。 设插值函数为1()L x ax b =+，由插值条件，建立线性方程组为 1 01 1a b a b e -?+=???+=? 解之得11 1a e b -?=-?=? 则11()(1)1L x e x -=-+ 因为(),()x x y x e y x e --'''=-= 所以，插值余项为 (1)(2) (2)011 ()()()()() (1)! 1()()2!1 ()()()2!1 (0)(1)((0,1))2n r x f x p x f x n f x f x x x x e x x ξξπξπξξ+-=-=+= =--=--∈ 所以 01 0101 ()max max (1) 2111248x r x e x x e ξξ-≤≤≤≤-≤-=??=。 2选用合适的三次插值多项式来近似计算f 和f 。 解：设三次插值多项式为230123()f x a a x a x a x =+++，由插值条件，建立方程组为 23012323 012323 01232301 23(0.1)(0.1)(0.1)0.9950.30.30.30.995 0.70.70.70.7651.1 1.1 1.10.454 a a a a a a a a a a a a a a a a ?+?-+?-+?-=?+?+?+?=??+?+?+?=??+?+?+?=?
即 012301230123 123012312301230.10.010.0010.9950.10.010.0010.9950.30.090.0270.9950.40.080.02800.70.490.3430.7650.80.480.344 1.761.1 1.21 1.3310.454a a a a a a a a a a a a a a a a a a a a a a a a a a -+-=-+-=??+++=++=??? +++=++=??+++=?12301231232330.40.720.9880.3110.10.010.0010.9950.40.080.02800.320.288 1.760.384 3.831a a a a a a a a a a a a a ??????++=-? -+-=??++=??? +=? ?-=-? 解之得 01 230.416.293.489.98 a a a a =??=-?? =-??=? 则所求的三次多项式为23()0.41 6.29 3.489.98f x x x x =--+。 所以 2323 (0.2)0.41 6.290.2 3.480.29.980.20.91 (0.8)0.41 6.290.8 3.480.89.980.8 1.74f f =-?-?+?=-=-?-?+?=- 3、设(0,1,2,,)i x i n =L 是 n+1个互异节点，证明： （1）0()(0,1,2,,)n k k i i i x l x x k n ===∑L ； （2）0 ()()0(0,1,2,,)n k i i i x x l x k n =-==∑L 。 证明： （1）由拉格朗日插值定理，以x 0,x 1,x 2,…x n 为插值节点，对y=f(x)=x k 作n 次插值，插值多项式为 0()()n n i i i p x l x y ==∑， 而y i =x i k , 所以0 ()()()n n k n i i i i i i p x l x y l x x ====∑∑ 同时，插值余项 (1)(1)11 ()()()()()()0(1)!(1)! n k n k n r x x p x f x x x n n ξξππ++=-= ==++ 所以0 ()n k k i i i l x x x ==∑ 结论得证。 （2）取函数()(),0,1,2,,k f x x t k n =-=L 对此函数取节点(0,1,2,,)i x i n =L ，则对应的插值多项式为
数值计算方法试题一 一、填空题（每空1分，共17分） 1、如果用二分法求方程在区间内的根精确到三位小数，需对分（）次。 2、迭代格式局部收敛的充分条件是取值在（）。 3、已知是三次样条函数，则 =( )，=（），=（）。 4、是以整数点为节点的Lagrange插值基函数，则 ( )，( )，当时( )。 5、设和节点则 和。 6、5个节点的牛顿-柯特斯求积公式的代数精度为，5个节点的求积公式最高代数精度为。 7、是区间上权函数的最高项系数为1的正交多项式族，其中，则。 8、给定方程组，为实数，当满足，且时，SOR迭代法收敛。 9、解初值问题的改进欧拉法是 阶方法。 10、设，当（）时，必有分解式，其中为下三角阵，当其对角线元素满足（）条件时，这种分解是唯一的。 二、二、选择题（每题2分） 1、解方程组的简单迭代格式收敛的充要条件是（）。（1）, (2) , (3) , (4) 2、在牛顿-柯特斯求积公式：中，当系数是负值时，公式的稳定性不能保证，所以实际应用中，当（）时的牛顿-柯特斯求积公式不使用。 （1），（2），（3），（4）， （1）二次；（2）三次；（3）四次；（4）五次 4、若用二阶中点公式求解初值问题，试问为保证该公式绝对稳定，步长的取值范围为（）。 (1), (2), (3), (4)
三、1、 2、（15 (1)(1) 试用余项估计其误差。 （2）用的复化梯形公式（或复化 Simpson公式）计算出该积分的近似值。 四、1、（15分）方程在附近有根，把方程写成三种不同的等价形式（1）对应迭代格式；(2)对应迭代格式；（3）对应迭代格式。判断迭代格式在的收敛性，选一种收敛格式计算附近的根，精确到小数点后第三位。选一种迭代格式建立Steffensen迭代法，并进行计算与前一种结果比较，说明是否有加速效果。 2、（8分）已知方程组，其中 ， （1）（1）列出Jacobi迭代法和Gauss-Seidel迭代法的分量形式。 （2）（2）求出Jacobi迭代矩阵的谱半径，写出SOR 迭代法。 五、1、（15分）取步长，求解初值问题用改进的欧拉法求的值；用经典的四阶龙格—库塔法求的值。 2、（8分）求一次数不高于4次的多项式使它满足 ,,,, 六、（下列2题任选一题，4分） 1、1、数值积分公式形如 （1）（1）试确定参数使公式代数精度尽量高；（2）设，推导余项公式，并估计误差。 2、2、用二步法 求解常微分方程的初值问题时，如何选择参数使方法阶数尽可能高，并求局部截断误差主项，此时该方法是几阶的。 数值计算方法试题二 一、判断题：（共16分，每小题２分） １、若是阶非奇异阵，则必存在单位下三角阵和上三角阵，使唯一成立。（）
1、解：将)(x V n 按最后一行展开，即知)(x V n 是n 次多项式。 由于 n i i i n n n n n i n x x x x x x x x x x V ...1...1... ......... ...... 1 )(21110 20 0---= ，.1,...,1,0-=n i 故知0)(=i n x V ，即110,...,,-n x x x 是)(x V n 的根。又)(x V n 的最高 次幂 n x 的系数为 )(...1...1... ...... .........1),...,,(101 1 21 11 2 2221 02001101j n i j i n n n n n n n n n n n x x x x x x x x x x x x x x V -== ∏-≤<≤-----------。 故知).)...()()(,...,,()(1101101------=n n n n x x x x x x x x x V x V 6、解：（1）设 .)(k x x f =当n k ,...,1,0=时，有.0)()1(=+x f n 对 )(x f 构造Lagrange 插值多项式， ),()(0 x l x x L j n j k j n ∑== 其 0)()! 1() ()()()(1)1(=+=-=++x w n f x L x F x R n n n n ξ， ξ介于j x 之间，.,...,1,0n j = 故 ),()(x L x f n =即 .,...,1,0,)(0 n k x x l x k j n j k j ==∑= 特别地，当0=k 时， 10) (=∑=n j x j l 。 （2） 0)()1(1) ()1()()(0000=-=??? ? ??-??? ? ??-=--=-===∑∑∑∑k j j i j i k j k i i j i i k j n j k i i j k n j j x x x x i k x l x x i k x l x x ）利用（。 7、证明：以b a ,为节点进行线性插值，得 )()()(1 b f a b a x a f b a b x x P --+--= 因 0)()(==b f a f ，故0)(1=x P 。而 ))()(("2 1 )()(1b x a x f x P x f --= -ξ，b a <<ξ。 故)("max )(8 122)("max )(max 2 2 x f a b a b x f x f b x a b x a b x a ≤≤≤≤≤≤-=??? ??-≤。 14、解：设 ))...()(()(21n n x x x x x x a x f ---=， k x x g =)(，记)() (1 ∏=-=n j j n x x x w ，则 ),()(x w a x f n n =).()(' j n n j x w a x f = 由差商的性质知 [])! 1()(1,..,,1) (' 1 )(')('1 211 11 -== ==-===∑∑∑ n g a x x x g a x w x a x w a x x f x n n n n n j j n k j n n j j n n k j n j j k j ξ， ξ介于n x x ,...,1之间。 当20-≤≤ n k 时，0)()1(=-ξn g ， 当 1-=n k 时，)!1()(1-=-n g n ξ， 故 ???-=-≤≤=-= --=∑1,,20,0)!1()(1) ('1 11 n k a n k n g a x f x n n n n j j k j ξ 16、解：根据差商与微商的关系，有 [] 1! 7! 7!7)(2,...,2,2)7(7 10===ξf f ， [ ] 0! 80 !8)(2,...,2,2)8(8 1 ===ξf f 。 （ 13)(47+++=x x x x f 是7次多项式， 故 ,!7)()7(=x f 0)()8(=x f ）。 25、解：（1） 右边= [][]dx x S x f x S dx x S x f b a b a ??-+-)(")(")("2)(")("2 = [] d x x S x f x S x S x S x f x f b a ?-++-)("2)(")("2)(")(")("2)(" 222 = [] d x x S x f b a ?-)(")(" 22 = [][]dx x S dx x f b a b a 2 2 )(")("??- =左边。 （2）左边= ? -b a dx x S x f x S ))(")(")(("
有限元分析概念 有限元法：把求解区域看作由许多小的在节点处相互连接的单元（子域）所构成，其模型给出基本方程的分片（子域）近似解，由于单元（子域）可以被分割成各种形状和大小不同的尺寸，所以它能很好地适应复杂的几何形状、复杂的材料特性和复杂的边界条件 有限元模型：它是真实系统理想化的数学抽象。由一些简单形状的单元组成，单元之间通过节点连接，并承受一定载荷。 有限元分析：是利用数学近似的方法对真实物理系统（几何和载荷工况）进行模拟。并利用简单而又相互作用的元素，即单元，就可以用有限数量的未知量去逼近无限未知量的真实系统。 线弹性有限元是以理想弹性体为研究对象的，所考虑的变形建立在小变形假设的基础上。在这类问题中，材料的应力与应变呈线性关系，满足广义胡克定律；应力与应变也是线性关系，线弹性问题可归结为求解线性方程问题，所以只需要较少的计算时间。如果采用高效的代数方程组求解方法，也有助于降低有限元分析的时间。 线弹性有限元一般包括线弹性静力学分析与线弹性动力学分析两方面。 非线性问题与线弹性问题的区别： 1）非线性问题的方程是非线性的，一般需要迭代求解； 2）非线性问题不能采用叠加原理； 3）非线性问题不总有一致解，有时甚至没有解。 有限元求解非线性问题可分为以下三类：
1）材料非线性问题 材料的应力和应变是非线性的，但应力与应变却很微小，此时应变与位移呈线性关系，这类问题属于材料的非线性问题。由于从理论上还不能提供能普遍接受的本构关系，所以，一般材料的应力与应变之间的非线性关系要基于试验数据，有时非线性材料特性可用数学模型进行模拟，尽管这些模型总有他们的局限性。在工程实际中较为重要的材料非线性问题有：非线性弹性（包括分段线弹性）、弹塑性、粘塑性及蠕变等。 2）几何非线性问题 几何非线性问题是由于位移之间存在非线性关系引起的。 当物体的位移较大时，应变与位移的关系是非线性关系。研究这类问题一般都是假定材料的应力和应变呈线性关系。它包括大位移大应变及大位移小应变问题。如结构的弹性屈曲问题属于大位移小应变问题，橡胶部件形成过程为大应变问题。 3）非线性边界问题 在加工、密封、撞击等问题中，接触和摩擦的作用不可忽视，接触边界属于高度非线性边界。 平时遇到的一些接触问题，如齿轮传动、冲压成型、轧制成型、橡胶减振器、紧配合装配等，当一个结构与另一个结构或外部边界相接触时通常要考虑非线性边界条件。 实际的非线性可能同时出现上述两种或三种非线性问题。
数值计算方法习题一（2） 习题二（6） 习题三（15） 习题四（29） 习题五（37） 习题六（62） 习题七（70） 2009．9，9
习题一 1．设x >0相对误差为2%4x 的相对误差。 解：由自变量的误差对函数值引起误差的公式： (())(())'()()()() f x x f x f x x f x f x δδ?= ≈得 （1）()f x = 11 ()()*2%1% 22x x δδδ≈ ===； （2）4 ()f x x =时 44 4 ()()'()4()4*2%8%x x x x x x δδδ≈ === 2．设下面各数都是经过四舍五入得到的近似数，即误差不超过最后一位的半个单位，试指出他们各有几位有效数字。 （1）12.1x =；（2）12.10x =；（3）12.100x =。 解：由教材9P 关于1212.m n x a a a bb b =±型数的有效数字的结论，易得上面三个数的有效 数字位数分别为：3，4，5 3．用十进制四位浮点数计算 （1）31.97+2.456+0.1352； （2）31.97+（2.456+0.1352） 哪个较精确？ 解：（1）31.97+2.456+0.1352 ≈2 1 ((0.3197100.245610)0.1352)fl fl ?+?+ =2 (0.3443100.1352)fl ?+ =0.3457210? （2）31.97+（2.456+0.1352） 2 1 (0.319710(0.245610))fl fl ≈?+? = 21 (0.3197100.259110)fl ?+? =0.34562 10? 易见31.97+2.456+0.1352=0.3456122 10?，故（2）的计算结果较精确。 4．计算正方形面积时，若要求面积的允许相对误差为1%，测量边长所允许的相对误差限为多少？
第一章绪论 习题一 1.设x>0,x*的相对误差为δ，求f(x)=ln x的误差限。解：求lnx的误差极限就是求f(x)=lnx的误差限，由公式(1. 2.4)有 已知x*的相对误差满足，而 ，故 即 2.下列各数都是经过四舍五入得到的近似值，试指出它们有几位有效数字，并给出其误差限与相对误差限。 解：直接根据定义和式(1.2.2)(1.2.3)则得 有5位有效数字，其误差限，相对误差限 有2位有效数字， 有5位有效数字， 3.下列公式如何才比较准确？ （1） （2）
解：要使计算较准确，主要是避免两相近数相减，故应变换所给公式。 （1） （2） 4.近似数x*=0.0310,是 3 位有数数字。 5.计算取，利用：式计算误差最小。 四个选项： 第二、三章插值与函数逼近 习题二、三 1. 给定的数值表 用线性插值与二次插值计算ln0.54的近似值并估计误差限. 解：仍可使用n=1及n=2的Lagrange插值或Newton插值,并应用误差估计（5.8）。线性插值时，用0.5及0.6两点，用Newton插值 误差限，因
，故 二次插值时，用0.5，0.6，0.7三点，作二次Newton插值 误差限 ，故 2. 在-4≤x≤4上给出的等距节点函数表，若用二次插值法求的近似值，要使误差不超过，函数表的步长h 应取多少? 解：用误差估计式（5.8）， 令 因 得 3. 若，求和.
解：由均差与导数关系 于是 4. 若互异，求 的值，这里p≤n+1. 解：，由均差对称性 可知当有 而当P＝n＋1时 于是得 5. 求证. 解：解：只要按差分定义直接展开得 6. 已知的函数表
习 题 三 解 答 1、用高斯消元法解下列方程组。 （1）1231231 22314254 27x x x x x x x x -+=?? ++=??+=?①②③ 解：?4②＋(-)①2，1 2 ?③＋(-)①消去第二、三个方程的1x ，得： 1232323231425313222 x x x x x x x ? ?-+=? -=???-=?④⑤⑥ 再由5 2)4 ?⑥＋(-⑤消去此方程组的第三个方程的2x ，得到三角方程组： 1232332314272184x x x x x x ? ?-+=? -=???-= ? 回代，得： 36x =-，21x =-，19x = 所以方程组的解为 (9,1,6)T x =-- 注意： ①算法要求，不能化简。化简则不是严格意义上的消元法，在算法设计上就多出了步骤。实际上，由于数值计算时用小数进行的，化简既是不必要的也是不能实现的。无论是顺序消元法还是选主元素消元法都是这样。 ②消元法要求采用一般形式,或者说是分量形式，不能用矩阵，以展示消元过程。 要通过练习熟悉消元的过程而不是矩阵变换的技术。 矩阵形式错一点就是全错，也不利于检查。 一般形式或分量形式： 1231231 22314254 27x x x x x x x x -+=?? ++=??+=?①②③ 矩阵形式 123213142541207x x x -?????? ??? ?= ??? ? ??? ???????
向量形式 123213142541207x x x -???????? ? ? ? ?++= ? ? ? ? ? ? ? ????????? ③必须是方程组到方程组的变形。三元方程组的消元过程要有三个方程组，不能变形出单一的方程。 ④消元顺序12x x →→L ，不能颠倒。按为支援在方程组中的排列顺序消元也是存储算法的要求。实际上，不按顺序消元是不规范的选主元素。 ⑤不能化简方程，否则系数矩阵会变化，也不利于算法设计。 （2）1231231231132323110 221x x x x x x x x x --=?? -++=??++=-? ①②③ 解：?23②＋( )①11，1 11 ?③＋(-)①消去第二、三个方程的1x ，得： 123232311323523569111111252414111111x x x x x x x ? --=?? ? -=? ? ? +=-??④⑤⑥ 再由25 11)5211 ?⑥＋(-⑤消去此方程组的第三个方程的2x ，得到三角方程组： 123233113235235691111111932235252x x x x x x ? ?--=? ? -=?? ? =-?? 回代，得： 32122310641 ,,193193193 x x x =- ==， 所以方程组的解为 41106223(,,)193193193T x =- 2、将矩阵 1020011120110011A ?? ? ?= ?- ???
《计算方法》期中复习试题 一、填空题： 1、已知3.1)3(,2.1)2(,0.1)1(===f f f ，则用辛普生（辛卜生）公式计算求得 ?≈3 1 _________ )(dx x f ,用三点式求得≈')1(f 。 答案：2.367，0.25 2、1)3(,2)2(,1)1(==-=f f f ，则过这三点的二次插值多项式中2 x 的系数为 ，拉 格朗日插值多项式为 。 答案：-1， )2)(1(21 )3)(1(2)3)(2(21)(2--------= x x x x x x x L 3、近似值*0.231x =关于真值229.0=x 有( 2 )位有效数字； 4、设)(x f 可微,求方程)(x f x =的牛顿迭代格式是( )； 答案 )(1)(1n n n n n x f x f x x x '--- =+ 5、对1)(3 ++=x x x f ,差商=]3,2,1,0[f ( 1 ),=]4,3,2,1,0[f ( 0 )； 6、计算方法主要研究( 截断 )误差和( 舍入 )误差； 7、用二分法求非线性方程 f (x )=0在区间(a ,b )内的根时，二分n 次后的误差限为 ( 1 2+-n a b )； 8、已知f (1)＝2，f (2)＝3，f (4)＝5.9，则二次Newton 插值多项式中x 2系数为( 0.15 )； 11、 两点式高斯型求积公式?1 d )(x x f ≈( ?++-≈1 )] 321 3()3213([21d )(f f x x f )，代数精度 为( 5 )； 12、 为了使计算 32)1(6 )1(41310-- -+-+ =x x x y 的乘除法次数尽量地少，应将该表达 式改写为 11 ,))64(3(10-= -++=x t t t t y ，为了减少舍入误差，应将表达式1999 2001-
有限单元法 有限元方法的基础是变分原理和加权余量法，其基本求解思想是把计算域划分为有限个互不重叠的单元，在每个单元内，选择一些合适的节点作为求解函数的插值点，将微分方程中的变量改写成由各变量或其导数的节点值与所选用的插值函数组成的线性表达式，借助于变分原理或加权余量法，将微分方程离散求解。采用不同的权函数和插值函数形式，便构成不同的有限元方法。有限元方法最早应用于结构力学，后来随着计算机的发展慢慢用于流体力学的数值模拟。在有限元方法中，把计算域离散剖分为有限个互不重叠且相互连接的单元，在每个单元内选择基函数，用单元基函数的线形组合来逼近单元中的真解，整个计算域上总体的基函数可以看为由每个单元基函数组成的，则整个计算域内的解可以看作是由所有单元上的近似解构成。在河道数值模拟中，常见的有限元计算方法是由变分法和加权余量法发展而来的里兹法和伽辽金法、最小二乘法等。根据所采用的权函数和插值函数的不同，有限元方法也分为多种计算格式。从权函数的选择来说，有配置法、矩量法、最小二乘法和伽辽金法，从计算单元网格的形状来划分，有三角形网格、四边形网格和多边形网格，从插值函数的精度来划分，又分为线性插值函数和高次插值函数等。不同的组合同样构成不同的有限元计算格式。对于权函数，伽辽金(Galerkin)法是将权函数取为逼近函数中的基函数；最小二乘法是令权函数等于余量本身，而内积的极小值则为对代求系数的平方误差最小；在配置法中，先在计算域内选取N个配置点。令近似解在选定的N个配置点上严格满足微分方程，即在配置点上令方程余量为0。插值函数一般由不同次幂的多项式组成，但也有采用三角函数或指数函数组成的乘积表示，但最常用的多项式插值函数。有限元插值函数分为两大类，一类只要求插值多项式本身在插值点取已知值，称为拉格朗日(Lagrange)多项式插值；另一种不仅要求插值多项式本身，还要求它的导数值在插值点取已知值，称为哈密特(Hermite)多项式插值。单元坐标有笛卡尔直角坐标系和无因次自然坐标，有对称和不对称等。常采用的无因次坐标是一种局部坐标系，它的定义取决于单元的几何形状，一维看作长度比，二维看作面积比，三维看作体积比。在二维有限元中，三角形单元应用的最早，近来四边形等参元的应用也越来越广。对于二维三角形和四边形电源单元，常采用的插值函数为有Lagrange插值直角坐标系中的线性插值函数及二阶或更高阶插值函数、面积坐标系中的线性插值函数、二阶或更高阶插值函数等。 对于有限元方法，其基本思路和解题步骤可归纳为 (1)建立积分方程，根据变分原理或方程余量与权函数正交化原理，建立与微分方程初边值问题等价的积分表达式，这是有限元法的出发点。 (2)区域单元剖分，根据求解区域的形状及实际问题的物理特点，将区域剖分为若干相互连接、不重叠的单元。区域单元划分是采用有限元方法的前期准备工作，这部分工作量比较大，除了给计算单元和节点进行编号和确定相互之间的关系之外，还要表示节点的位置坐标，同时还需要列出自然边界和本质边界的节点序号和相应的边界值。 (3)确定单元基函数，根据单元中节点数目及对近似解精度的要求，选择满足一定插值条件的插值函数作为单元基函数。有限元方法中的基函数是在单元中选取的，由于各单元具有规则的几何形状，在选取基函数时可遵循一定的法则。 (4)单元分析：将各个单元中的求解函数用单元基函数的线性组合表达式进行逼近；再将近似函数代入积分方程，并对单元区域进行积分，可获得含有待定系数(即单元中各节点的参数值)的代数方程组，称为单元有限元方程。 (5)总体合成：在得出单元有限元方程之后，将区域中所有单元有限元方程按一定法则进行累加，形成总体有限元方程。
1.什么是边缘计算？ 在IIoT的背景下，“边缘”是指靠近数据源的计算基础设施，例如工业机器（例如风力涡轮机，磁共振（MR）扫描仪，海底防喷器）），工业控制器如SCADA系统和时间序列数据库汇总来自各种设备和传感器的数据。这些设备通常远离云中可用的集中式计算。 边缘计算是在靠近物或数据源头的网络边缘侧，融合网络、计算、存储、应用核心能力的开放平台。边缘计算与云计算互相协同，共同助力各行各业的数字化转型。它就近提供智能互联服务，满足行业在数字化变革过程中对实时业务、业务智能、数据聚合与互操作、安全与隐私保护等方面的关键需求。 到目前为止，边缘计算的作用主要用于摄取，存储，过滤和发送数据到云系统。然而，我们正处于一个时间点，这些计算系统正在包装更多的计算，存储和分析功能，以消耗并对机器位置的数据采取行动。这种能力对于工业组织来说将是非常有价值的 - 这是不可或缺的。 2.这对工业带来的价值 行业权威人士已经计算出，数以千计的连接事物会从不同的来源产生大量的数据。根据国际电信联盟电信标准分局ITU-T的研究报告，到2020年，每个人每秒将产生的数据，IoT可穿戴设备的出货量将达到亿。IDC也发布了相关预测，到2018年，50%的物联网网络将面临网络带宽的限制，40%的数据需要在网络边缘侧分析、处理与储存，到2025年，这一数字将超过50%。管理咨询公司麦肯锡公司估计，到2025年，工业物联网（IIoT）将创造价值万亿的市场规模。工业物联网将思想和机器结合在一起，将人们与加速数字产业转型的机器数据相结合。 通过将大数据，高级分析和机器学习应用于运营，工业可以减少计划外停机时间，提高资产性能，降低维护成本，并为从机床数据中获取未开发价值的新业务模式开拓潜力。 过去几年来，工业组织已经开始将云计算融入业务，从大量数据中获取洞察力，帮助实现关键业务成果，包括减少意外停机，提高生产效率，降低能耗等。云计算仍然通过工业物联网来实现新的性能水平发挥关键作用，因为它需要大量的计算能力来有效地管理来自机器的庞大数据量。 但是随着更多的计算，存储和分析能力被捆绑到更靠近数据源的较小设备中，即工业机器 - 边缘计算将有助于边缘处理实现工业物联网的承诺。 虽然这个概念不是新的，但是有几个关键的驱动力使它成为今天更可行的现实：·计算和传感器的成本继续下滑， ·在较小尺寸的设备（如网关或传感器集线器）中执行的更多计算能力， ·来自机器和/或环境的日益增长的数据（例如天气或市场定价）， ·现代机器学习与分析。 这些因素有助于公司将大量数据转化为具有洞察力和智慧的行动。 对于工业组织来说，这种技术在以下用例中将变得至关重要： ·低/间歇连接（如远程位置） o将数据传输到云的带宽和相关的高成本 o低延迟，例如机器洞察和启动之间的闭环相互作用（即在机器上采取动作）
数值分析作业答案 插值法 1、当x=1，-1，2时，f(x)=0，-3，4，求f(x)的二次插值多项式。 （1）用单项式基底。 （2）用Lagrange插值基底。 （3）用Newton基底。 证明三种方法得到的多项式是相同的。 解：（1）用单项式基底 设多项式为: ， 所以： 所以f(x)的二次插值多项式为： （2）用Lagrange插值基底 Lagrange插值多项式为： 所以f(x)的二次插值多项式为： (3) 用Newton基底: 均差表如下： xk f(xk) 一阶均差二阶均差 1 0 -1 -3 3/2 2 4 7/ 3 5/6 Newton插值多项式为： 所以f(x)的二次插值多项式为： 由以上计算可知，三种方法得到的多项式是相同的。 6、在上给出的等距节点函数表，若用二次插值求ex的近似值，要使截断误差不超过10-6，问使用函数表的步长h应取多少？ 解：以xi-1,xi,xi+1为插值节点多项式的截断误差，则有 式中 令得 插值点个数
是奇数，故实际可采用的函数值表步长 8、，求及。 解：由均差的性质可知，均差与导数有如下关系： 所以有： 15、证明两点三次Hermite插值余项是 并由此求出分段三次Hermite插值的误差限。 证明：利用[xk，xk+1]上两点三次Hermite插值条件 知有二重零点xk和k+1。设 确定函数k(x): 当或xk+1时k(x)取任何有限值均可； 当时，，构造关于变量t的函数 显然有 在[xk,x][x,xk+1]上对g(x)使用Rolle定理，存在及使得 在，，上对使用Rolle定理，存在，和使得 再依次对和使用Rolle定理，知至少存在使得 而，将代入，得到 推导过程表明依赖于及x 综合以上过程有： 确定误差限： 记为f(x)在[a,b]上基于等距节点的分段三次Hermite插值函数。在区间[xk，xk+1]上有 而最值 进而得误差估计： 16、求一个次数不高于4次的多项式，使它满足，，。
比较详细的数值分析课后习题答案
0.1算法 1、 （p.11，题1）用二分法求方程013 =--x x 在[1,2]的近似根，要求误差不超过 10-3. 【解】 由二分法的误差估计式31 1*102 1 2||-++=≤=-≤ -εk k k a b x x ，得到100021≥+k .两端取自然对数得96.812ln 10 ln 3≈-≥ k ，因此取9=k ，即至少需 2、（p.11，题2） 证明方程210)(-+=x e x f x 在区间[0,1]有唯一个实根；使用二 分法求这一实根，要求误差不超过2102 1 -?。 【解】 由于210)(-+=x e x f x ，则)(x f 在区间[0,1]上连续，且 012010)0(0<-=-?+=e f ，082110)1(1>+=-?+=e e f ，即0)1()0(+=x e x f ，即)(x f 在区间[0,1]上是单调的，故)(x f 在区间[0,1]有唯一实根.
由二分法的误差估计式21 1*1021 2 12||-++?=≤=-≤-εk k k a b x x ，得到1002≥k .两端取自然对数得6438.63219.322 ln 10 ln 2=?≈≥ k ，因此取7=k ，即至少需二分 0.2误差 1．（p.12，题8）已知e=2.71828…，试问其近似值7.21=x ，71.22=x ，x 2=2.71， 718.23=x 各有几位有效数字？并给出它们的相对误差限。 【解】有效数字： 因为111021 05.001828.0||-?= <=- x e ，所以7.21=x 有两位有效数字； 因为1 2102105.000828.0||-?=<=- x e ，所以71.22=x 亦有两位有效数字； 因为3 3102 10005.000028.0||-?=<=- x e ，所以718.23=x 有四位有效数字； %85.17.205 .0||111=<-= x x e r ε； %85.171 .205 .0||222=<-= x x e r ε；
1.1 有限单元法中“离散”的含义是什么？有限单元法是如何将具有无限自由度的连续介质问题转变成有限自由度问题的？位移有限元法的标准化程式是怎样的？ （1）离散的含义即将结构离散化，即用假想的线或面将连续体分割成数目有限的单元，并在其上设定有限个节点；用这些单元组成的单元集合体代替原来的连续体，而场函数的节点值将成为问题的基本未知量。 （2）给每个单元选择合适的位移函数或称位移模式来近似地表示单元内位移分布规律，即通过插值以单元节点位移表示单元内任意点的位移。因节点位移个数是有限的，故无限自由度问题被转变成了有限自由度问题。 （3）有限元法的标准化程式：结构或区域离散，单元分析，整体分析，数值求解。 1.3 单元刚度矩阵和整体刚度矩阵各有哪些性质？各自的物理意义是什么?两者有何区别?单元刚度矩阵的性质：对称性、奇异性（单元刚度矩阵的行列式为零）。整体刚度矩阵的性质：对称性、奇异性、稀疏性。单元 Kij 物理意义 Kij 即单元节点位移向量中第 j 个自由度发生单位位移而其他位移分量为零时，在第 j 个自由度方向引起的节点力。整体刚度矩阵 K 中每一列元素的物理意义是：要迫使结构的某节点位移自由度发生单位位移，而其他节点位移都保持为零的变形状态，在所有个节点上需要施加的节点荷载。 2.2 什么叫应变能？什么叫外力势能？试叙述势能变分原理和最小势能原理，并回答下述问题：势能变分原理代表什么控制方程和边界条件？其中附加了哪些条件？ (1)在外力作用下，物体内部将产生应力σ和应变ε，外力所做的功将以变形能的形式储存起来，这种能量称为应变能。 (2)外力势能就是外力功的负值。 (3)势能变分原理可叙述如下：在所有满足边界条件的协调位移中，那些满足静力平衡条件的位移使物体势能泛函取驻值，即势能的变分为零 δ∏p=δ Uε+δV=0 此即变分方程。对于线性弹性体，势能取最小值，即 δ2∏P=δ2Uε+δ2V≥0 此时的势能变分原理就是著名的最小势能原理。 势能变分原理代表平衡方程、本构方程和应力边界条件，其中附加了几何方程和位移边界条件。 2.3 什么是强形式？什么是弱形式？两者有何区别？建立弱形式的关键步骤是什么？ 等效积分形式通过分部积分，称式 ∫ΩCT(v)D(u)dΩ+∫ΓET(v)F(u)dΓ 为微分方程的弱形式，相对而言，定解问题的微分方程称为强形式。 区别：弱形式得不到解析解。建立弱形式的关键步骤：对场函数要求较低阶的连续性。2.4 为了使计算结果能够收敛于精确解，位移函数需要满足哪些条件？为什么？ 只要位移函数满足两个基本要求，即完备性和协调性，计算结果便收敛于精确解。 2.6 为什么采用变分法求解通常只能得到近似解？变分法的应用常遇到什么困难？Ritz 法收敛的条件是什么？ （1）在 Ritz 法中，N 决定了试探函数的基本形态，待定参数使得场函数具有一定的任意性。如果真实场函数包含在试探函数之内，则变分法得到的解答是精确的；如果试探函数取自完全的函数序列，则当项数不断增加时，近似解将趋近于精确解。然而，通常情况下试探函数不会将真实场函数完全包含在内，实际计算时也不可能取无穷多项。因此，试探函数只能是真实场函数的近似。可见，变分法就是在某个假定的范围内找出最佳解答，近似性就源于此。 （2）采用变分法近似求解，要求在整个求解区域内预先给出满足边界条件的场函数。通常情况下这是不可能的，因而变分法的应用受到了限制。 （3）Ritz 法的收敛条件是要求试探函数具有完备性和连续性，也就是说，如果试探函数满足完备性和连续性的要求，当试探函数的项数趋近于无穷时，则 Ritz 法的近似解将趋近于数学微分方程的精确解。 3.1 构造单元形函数有哪些基本原则？ 形函数是定义于单元内坐标的连续函数。单元位移函数通常采用多项式，其中的待定常数应该与单元节点自由度数相等。为满足完备性要求，位移函数中必须包括常函数和一次式，即完全一次多项式。多项式的选取应由低阶到高阶，尽量选择完全多项式以提高单元的精度。若由于项数限制而不能选取完全多项式时，也应使完全多项式具有坐标的对称性，并且一
《数值计算方法》复习试题 一、填空题： 1、 ?? ??? ?????----=410141014A ，则A 的LU 分解为 A ? ???????? ???=????????? ?? ?。 答案： ?? ????????--??????????--=1556141501 4115401411A 2、已知3.1)3(,2.1)2(, 0.1)1(===f f f ，则用辛普生（辛卜生）公式计算求 得?≈3 1 _________ )(dx x f ,用三点式求得≈')1(f 。 答案：2.367，0.25 3、1)3(,2)2(, 1)1(==-=f f f ，则过这三点的二次插值多项式中2x 的系数 为 ，拉格朗日插值多项式为 。 答案：-1， )2)(1(21 )3)(1(2)3)(2(21)(2--------= x x x x x x x L 4、近似值*0.231x =关于真值229.0=x 有( 2 )位有效数字； 5、设)(x f 可微,求方程)(x f x =的牛顿迭代格式是( )； 答案 )(1)(1n n n n n x f x f x x x '--- =+ 6、对 1)(3++=x x x f ,差商=]3,2,1,0[f ( 1 ),=]4,3,2,1,0[f ( 0 )； 7、计算方法主要研究( 截断 )误差和( 舍入 )误差； 8、用二分法求非线性方程f (x )=0在区间(a ,b )内的根时，二分n 次后的误差限为( 1 2+-n a b )； 9、求解一阶常微分方程初值问题y '= f (x ,y )，y (x 0)=y 0的改进的欧拉公