有限元法的预备知识.ppt

第二章有限元法的预备知识 矩阵 线性代数和MATLAB 在数学上 有限元法通过变分法或加权残数法把微分方程化为近似的代数方程组进行求解 如果采用矩阵形式进行中间过程的推导 就能使有限元公式简洁紧凑 为此 我们必须掌握矩阵和线性代数的基本知识 了解矩阵的详细操作 这非常有助于阅读和理解有限元法的基本理论并组织有效的算法 在叙述矩阵 线性代数基本知识的同时 将介绍使用MATLAB实现矩阵定义 运算等基本操作 这有助于增加感性认识 加深理解 矩阵简介 A 1 40 0 1 0 0 20 01 4 0 40 4 1 0 0 41 80 6 0 20 40 62 4 A 1 40000 1 0000 0 200001 4000 0 40000 4000 1 0000 0 40001 80000 6000 0 20000 40000 60002 4000 在MATLAB的命令窗口中输入 MATLAB在命令窗口显示如下 表示刚才的输入正确 已经被MATLAB接受 下次可以引用矩阵A 元素之间可以用空格或者逗号隔开 行与行之间可以用回车或者分号隔开 列向量b可以这样输入 b 1 0 0 0 0 0 0 0 特殊矩阵 1 零矩阵 MATLAB中定义零矩阵十分方便 比如定义一个4 3的零矩阵的命令是 A zeros 4 3 A 000000000000 MATLAB在命令窗口显示如下 命令zeros接受可接受任意个参数 一般我们使用的都是2维的矩阵 前一个参数表示行数 后一个是列数 2 对称矩阵 目前 在MATLAB中没有定义对称矩阵的函数 或者说没有对称矩阵概念 这是MATLAB中一个比较大的缺憾 MATLAB的函数transpose是专门对矩阵进行转置操作的 例如 A 1234 31 A 1234 31 transpose A ans 142 331 函数transpose可以用一个对应的操作符单引号 来代替 例如 A ans 142 331 这里需要指出的是 函数transpose与单引号对于实矩阵的操作是完全一样的 而对于复矩阵则不同 因为复矩阵在数学上的转置定义不仅包括行列互换 还要执行共轭操作 函数transpose只执行行列互换的操作 而单引号则是数学意义上的转置 3 单位矩阵 MATLAB中用eye来定义单位矩阵 下面的命令定义了一个三阶的单位矩阵 eye 3 ans 100010001 函数eye也能接受多个参数 当它们不相等时 给出的是对角线元素为1其余元素为零的非方矩阵 例如 eye 5 4 ans 10000100001000010000 为什么不用I 4 带状稀疏矩阵 对称带状矩阵的一维存储 在MATLAB中 我们可以用命令sparse定义一个稀疏矩阵 比如 A sparse 4000 4000 命令窗口中显示A Allzerosparse 4000 by 4000用来表示A是一个4000 4000的元素全部为零的稀疏矩阵 对矩阵元素的引用跟普通的矩阵一样 矩阵的运算 1 相等 2 加法与减法 MATLAB中矩阵的运算跟单个数值的运算具有相同的形式 比如 A 123 4 31 B 11 1 32 5 C 3 14 263 A Bans 2321 1 4 B Aans 2321 1 4 A B Cans 126 15 1 A B C ans 126 15 1 操作方便 数乘矩阵 A 123 4 31 B 11 1 263 k 4 l 2 k Aans 481216 124 A kans 481216 124 k l A ans 8162432 248 k l Aans 8162432 248 k l Aans 6121824 186 k A l Aans 6121824 186 k A B ans 812881216 k A k Bans 812881216 操作方便 矩阵相乘 矩阵相乘的性质 A 1 2 B 34 A Bans 3468 B Aans 11 我们用MATLAB来验证上式 即 矩阵相乘的性质 下面用MATLAB验证矩阵乘法的分配律 先定义三个随机矩阵 A rand 4 3 B rand 4 3 C rand 3 2 然后验证乘法的分配律 即 A B Cans 0 80511 80230 63361 20090 56840 99280 46810 9461 A C B Cans 0 80511 80230 63361 20090 56840 99280 46810 9461 操作方便 矩阵不能约去 但是 矩阵积的转置 如果 那么 逆矩阵 如果 那么 矩阵的求逆对于一般的程序语言必须编写一段复杂的程序 根据矩阵的特点 选择合适的算法 十分麻烦 而对于MATLAB 则是轻而易举 一个inv函数就足够了 例如 A rand 4 4 A 0 01530 46600 20260 68130 74680 41860 67210 37950 44510 84620 83810 83180 93180 52520 01960 5028 B inv A B 0 22031 3026 1 11100 5563 5 0783 5 45525 56051 79910 15781 6202 0 1151 1 24614 88953 2202 3 7443 0 8724上面的矩阵B就是A的逆矩阵 我们可以通过下面的命令可以验证它们的互逆关系 A Bans 1 000000001 000000001 000000001 0000A和B相乘的结果是单位阵 因此它们是互逆的 矩阵分块 MATLAB中分块矩阵是通过冒号来实现的 假设A是一个5 5的矩阵 那么A 2 表示矩阵A的第二行A 3 表示矩阵A的第三列A 2 4 表示由矩阵A的第2 3和4行元素组成的一个3 5的矩阵 A 3 5 表示由矩阵A的第3 4和5列元素组成的一个5 3的矩阵A 2 4 3 5 表示由矩阵A的第2行到第4行 第3列到第5列元素组成的一个3 3的矩阵 A 1 2 5 表示由矩阵A的第1 3和5行元素组成的一个3 5的矩阵 A rand 5 5 A 0 70950 68220 37840 89980 34200 42890 30280 86000 82160 28970 30460 54170 85370 64490 34120 18970 15090 59360 81800 53410 19340 69790 49660 66020 7271 A 2 ans 0 42890 30280 86000 82160 2897 A 3 ans 0 37840 86000 85370 59360 4966 A 2 4 ans 0 42890 30280 86000 82160 28970 30460 54170 85370 64490 34120 18970 15090 59360 81800 5341 A 0 70950 68220 37840 89980 34200 42890 30280 86000 82160 28970 30460 54170 85370 64490 34120 18970 15090 59360 81800 53410 19340 69790 49660 66020 7271 A 3 5 ans 0 37840 89980 34200 86000 82160 28970 85370 64490 34120 59360 81800 53410 49660 66020 7271 A 2 4 3 5 ans 0 86000 82160 28970 85370 64490 34120 59360 81800 5341 A 1 2 5 ans 0 70950 68220 37840 89980 34200 30460 54170 85370 64490 34120 19340 69790 49660 66020 7271 矩阵的迹和行列式 MATLAB中用两个函数trace和det分别求取矩阵的迹和行列式值 A 346302 2500 230004 A 346302 2500 230004 trace A ans 7 det A ans 48 太简单了 MATLAB语言程序设计初步 ifI JA I J 2 elseifabs I J 1A I J 1 elseA I J 0 end IF语句 MATLAB与C C 十分相似 helpifIFstatementcondition ThegeneralformoftheIFstatementisIFexpressionstatementsELSEIFexpressionstatementsELSEstatementsENDThestatementsareexecutediftherealpartoftheexpressionhasallnon zeroelements TheELSEandELSEIFpartsareoptional ZeroormoreELSEIFpartscanbeusedaswellasnestedIF s Theexpressionisusuallyoftheformexprropexprwhereropis or SWITCH语句 method Bilinear switchlower method case linear bilinear disp Methodislinear case cubic disp Methodiscubic case nearest disp Methodisnearest otherwisedisp Unknownmethod end helpswitchSwitchamongseveralcasesbasedonexpression ThegeneralformoftheSWITCHstatementis SWITCHswitch exprCASEcase expr statement statementCASE case expr1 case expr2 case expr3 statement statement OTHERWISE statement statementENDThestatementsfollowingthefirstCASEwheretheswitch exprmatchesthecase exprareexecuted Whenthecaseexpressionisacellarray asinthesecondcaseabove thecase exprmatchesifanyoftheelementsofthecellarraymatchtheswitchexpression IfnoneofthecaseexpressionsmatchtheswitchexpressionthentheOTHERWISEcaseisexecuted ifitexists OnlyoneCASEisexecutedandexecutionresumeswiththestatementaftertheEND Theswitch exprcanbeascalarorastring Ascalarswitch exprmatchesacase exprifswitch expr case expr Astringswitch exprmatchesacase exprifstrcmp switch expr case expr returns1 true OnlythestatementsbetweenthematchingCASEandthenextCASE OTHERWISE orENDareexecuted UnlikeC theSWITCHstatementdoesnotfallthrough soBREAKsareunnecessary FOR语句 N 10 forI 1 NforJ 1 NA I J 1 I J 1 endend helpforFORRepeatstatementsaspecificnumberoftimes ThegeneralformofaFORstatementis FORvariable expr statement statementENDThecolumnsoftheexpressionarestoredoneatatimeinthevariableandthenthefollowingstatements uptotheEND areexecuted TheexpressionisoftenoftheformX Y inwhichcaseitscolumnsaresimplyscalars Someexamples assumeNhasalreadybeenassignedavalue WHILE语句 E 0 A F E eye size E N 1 whilenorm E F E 1 0E E F F A F N N N 1 end helpwhileWHILERepeatstatementsanindefinitenumberoftimes ThegeneralformofaWHILEstatementis WHILEexpressionstatementsENDThestatementsareexecutedwhiletherealpartoftheexpressionhasallnon zeroelements Theexpressionisusuallytheresultofexprropexprwhereropis or TheBREAKstatementcanbeusedtoterminatetheloopprematurely sinpn1 1 0 fork 2 15n 2 k sinpn sqrt 2 2 sqrt 1 sqrt 1 sinpn1 2 sinpn1 sinpn pn n sinpn fprintf n 7dpi 10f n n pn end 一个简单的例子 计算圆周率的近似方法如下 用MATLAB实现上述公式的代码为 FUNCTION语句 function mean stdev stat x mainfunctionn length x mean avg x n stdev sqrt sum x avg x n 2 n functionmean avg x n AVGsub functionmean sum x n 假设硬盘上有文件stat m 其中有这样两个函数 函数stat可以在其他m文件中被调用 而