当前位置:文档之家› 列主元三角分解法在matlab中的实现

列主元三角分解法在matlab中的实现

列主元三角分解法在matlab中的实现
列主元三角分解法在matlab中的实现

列主元三角分解法在matlab中的实现

列主元三角分解法在matlab中的实现

摘要:介绍了M atlab语言并给出用M atlab语言实现线性方程组的列主元三角分解法,其有效性已在计算机实现中得到了验证。

关键词:M atlab语言;高斯消去法;列主元三角分解法

0前言

M atlab是M atrix Laboratory(矩阵实验室)的缩写,它是由美国M athwork公司于1967年推出的软件包,现已发展成为一种功能强大的计算机语言。它编程简单,使用方便,在M a tlab环境下数组的操作与数的操作一样简单,进行数学运算可以像草稿纸一样随心所欲,使计算机兼备高级计算器的优点。M atlab语言具有强大的矩阵和向量的操作功能,是Fo rtran和C语言无法比拟的;M a tlab语言的函数库可任意

扩充;语句简单,内涵丰富;还具有二维和三维绘图功能且使用方便,特别适用于科学和工程计算。

在科学和工程计算中,应用最广泛的是求解线性方程组的解,一般可用高斯消去法求解,如果系数矩阵不满足高斯消去法在计算机上可行的条件,那么消元过程中可能会出现零主元或小主元,消元或不可行或数值不稳定,解决办法就是对方程组进行行交换或列交换来消除零主元或小主元,这就是选主元的思想。

1 定义

列主元三角分解:如果A为非奇异矩阵,则存在排列矩阵P,使PA=LU,其中L为单位下三角矩阵,U为上三角阵。列主元三角分角法

是对直接三角分解法的一种改进,主要目的和列主元高斯消元法一样,就是避免小数作为分母项.

2 算法概述

列主元三角分解法和普通三角分解法基本上类似,所不同的是在构造Gauss 变换前,先在对应列中选择绝对值最大的元素(称为列主元),然后实施初等行交换将该元素调整到矩阵对角线上。 例如第)1,,2,1(-=n k Λ步变换叙述如下:

选主元:确定p 使{}1)1(

max -≤≤-=k ik n

i k k pk a a ; 行交换:将矩阵的第k 行和第p 行上的元素互换位置,即

实施Gauss 变换:通过初行变换,将列主对角线以下的元素消为零.即

3 列主元三角分解在matlab 中的实现

例:对矩阵5

132523

21 A 进行LPU 分解

其程序如下:

function [l,u,p]=mylu(A)

[m,n]=size(A);

if m~=n

error('矩阵不是方阵')

return

end

if det(A)==0

error('矩阵不能被三角分解')

end

u=A;p=eye(m);l=eye(m);

for i=1:m

for j=i:m

t(j)=u(j,i);

for k=1:i-1

t(j)=t(j)-u(j,k)*u(k,i);

end

end

a=i;b=abs(t(i));

for j=i+1:m

if b

b=abs(t(j)); a=j;

end

end

if a~=i

for j=1:m

c=u(i,j);

u(i,j)=u(a,j); u(a,j)=c;

end

for j=1:m

c=p(i,j);

p(i,j)=p(a,j); p(a,j)=c;

end

c=t(a);

t(a)=t(i);

t(i)=c;

end

u(i,i)=t(i);

for j=i+1:m

u(j,i)=t(j)/t(i);

end

for j=i+1:m

for k=1:i-1

u(i,j)=u(i,j)-u(i,k)*u(k,j);

end

end

end

l=tril(u,-1)+eye(m);

u=triu(u,0);

控制命令为:

A=[1 2 3;2 5 2;3 1 5];

[l,u,p]=mylu(A)

结果为:

l =

1.0000 0 0

0.6667 1.0000 0

0.3333 0.3846 1.0000

u =

3.0000 1.0000 5.0000

0 4.3333 -1.3333

0 0 1.8462

p =

0 0 1

0 1 0

1 0 0

4 小结

在数值计算中,列主元三角分解法是求解线性方程组一个很重要的方法,而用MATLAB可以简单便捷的实现该算法,从而轻松得到

线性方程组的解。

参考文献:

[1]蒲俊等,吉家锋,伊良忠.MATLAB6.0数学手册[M].上海:浦东电子出版社,2002.35-40.

[2]萧树铁.数学实验[M].北京:高等教育出版社,1999.130-139.

[3]李庆扬,王能超,易大义.数值分析(第7、8章)[M].武汉:华中理工大学出版社,1987.

[4]李丽,王振领.MATLAB工程计算及应用[M].北京:人民邮电出版社,2001.169-172.

matlab有限元分析实例

MATLAB: MATLAB是美国MathWorks公司出品的商业数学软件,用于数据分析、无线通信、深度学习、图像处理与计算机视觉、信号处理、量化金融与风险管理、机器人,控制系统等领域。 MATLAB是matrix&laboratory两个词的组合,意为矩阵工厂(矩阵实验室),软件主要面对科学计算、可视化以及交互式程序设计的高科技计算环境。它将数值分析、矩阵计算、科学数据可视化以及非线性动态系统的建模和仿真等诸多强大功能集成在一个易于使用的视窗环境中,为科学研究、工程设计以及必须进行有效数值计算的众多科学领域提供了一种全面的解决方案,并在很大程度上摆脱了传统非交互式程序设计语言(如C、Fortran)的编辑模式。 MATLAB和Mathematica、Maple并称为三大数学软件。它在数学类科技应用软件中在数值计算方面首屈一指。MATLAB可以进行矩阵运算、绘制函数和数据、实现算法、创建用户界面、连接其他编程语言的程序等。MATLAB的基本数据单位是矩阵,它的指令表达式与数学、工程中常用的形式十分相似,故用MATLAB来解算问题要比用C,FORTRAN等语言完成相同的事情简捷得多,并且MATLAB也吸收了像Maple等软件的优点,使MATLAB成为一个强大的数学软件。在新的版本中也加入了对C,FORTRAN,C++,JAVA的支持。 MATLAB有限元分析与应用:

《MATLAB有限元分析与应用》是2004年4月清华大学出版社出版的图书,作者是卡坦,译者是韩来彬。 内容简介: 《MATLAB有限元分析与应用》特别强调对MATLAB的交互应用,书中的每个示例都以交互的方式求解,使读者很容易就能把MATLAB用于有限分析和应用。另外,《MATLAB有限元分析与应用》还提供了大量免费资源。 《MATLAB有限元分析与应用》采用当今在工程和工程教育方面非常流行的数学软件MATLAB来进行有限元的分析和应用。《MATLAB有限元分析与应用》由简单到复杂,循序渐进地介绍了各种有限元及其分析与应用方法。书中提供了大量取自机械工程、土木工程、航空航天工程和材料科学的示例和习题,具有很高的工程应用价值。

列主元高斯消去法和列主元三角分解法解线性方程

计算方法实验报告1 【课题名称】 用列主元高斯消去法和列主元三角分解法解线性方程 【目的和意义】 高斯消去法是一个古老的求解线性方程组的方法,但由它改进得到的选主元的高斯消去法则是目前计算机上常用的解低阶稠密矩阵方程组的有效方法。 用高斯消去法解线性方程组的基本思想时用矩阵行的初等变换将系数矩阵A 约化为具有简单形式的矩阵(上三角矩阵、单位矩阵等),而三角形方程组则可以直接回带求解 用高斯消去法解线性方程组b Ax =(其中A ∈Rn ×n )的计算量为:乘除法运算步骤为 32(1)(1)(21)(1)(1)262233n n n n n n n n n n n MD n ----+= +++=+-,加减运算步骤为 (1)(21)(1)(1)(1)(25) 6226 n n n n n n n n n n AS -----+= ++= 。相比之下,传统的克莱姆 法则则较为繁琐,如求解20阶线性方程组,克莱姆法则大约要19 510?次乘法,而用高斯消去法只需要3060次乘除法。 在高斯消去法运算的过程中,如果出现abs(A(i,i))等于零或过小的情况,则会导致矩阵元素数量级严重增长和舍入误差的扩散,使得最后的计算结果不可靠,所以目前计算机上常用的解低阶稠密矩阵方程的快速有效的方法时列主元高斯消去法,从而使计算结果更加精确。 2、列主元三角分解法 高斯消去法的消去过程,实质上是将A 分解为两个三角矩阵的乘积A=LU ,并求解Ly=b 的过程。回带过程就是求解上三角方程组Ux=y 。所以在实际的运算中,矩阵L 和U 可以直接计算出,而不需要任何中间步骤,从而在计算过程中将高斯消去法的步骤进行了进一步的简略,大大提高了运算速度,这就是三角分解法 采用选主元的方式与列主元高斯消去法一样,也是为了避免除数过小,从而保证了计算的精确度 【计算公式】 1、 列主元高斯消去法 设有线性方程组Ax=b ,其中设A 为非奇异矩阵。方程组的增广矩阵为 第1步(k=1):首先在A 的第一列中选取绝对值最大的元素 1l a ,作为第一步的主元素: 111211212222112[,]n n n l n nn n a a a a b a a a b a a a b ?? ???? ?? =?????? ?? ????a b

Matlab有限元分析操作基础

Matlab 有限元分析20140226 为了用Matlab 进行有限元分析,首先要学会Matlab 基本操作,还要学会使用Matlab 进行有限元分析的基本操作。 1. 复习:上节课分析了弹簧系统 x 推导了系统刚度矩阵 11221 21200k k k k k k k k -????-????--+??

2. Matlab有限元分析的基本操作 (1)单元划分(选择何种单元,分成多少个单元,标号)(2)构造单元刚度矩阵(列出…) (3)组装系统刚度矩阵(集成整体刚度矩阵) (4)引入边界条件(消除冗余方程) (5)解方程 (6)后处理(扩展计算)

3. Matlab有限元分析实战【实例1】

分析: 步骤一:单元划分

步骤二:构造单元刚度矩阵 >>k1=SpringElementStiffness(100) >>…?

步骤三:构造系统刚度矩阵 a) 分析SpringAssemble库函数function y = SpringAssemble(K,k,i,j) % This function assembles the element stiffness % matrix k of the spring with nodes i and j into the % global stiffness matrix K. % function returns the global stiffness matrix K % after the element stiffness matrix k is assembled. K(i,i) = K(i,i) + k(1,1); K(i,j) = K(i,j) + k(1,2); K(j,i) = K(j,i) + k(2,1); K(j,j) = K(j,j) + k(2,2); y = K; b) K是多大矩阵? 今天的系统刚度矩阵是什么? 因为 11 22 1212 k k k k k k k k - ?? ?? - ????--+ ?? 所以 1000100 0200200 100200300 - ?? ?? - ?? ?? -- ?? ?

2.4直接三角分解法

§4 直接三角分解法 一、教学设计 1.教学内容:Doolittle 分解法、Crout 分解法,紧凑格式的Doolittle 分解法、部分选主元的Doolittle 分解法。 2.重点难点:紧凑格式的Doolittle 分解法、部分选主元的Doolittle 分解法。 3.教学目标:了解直接三角分解法的基本思想,掌握基本三角分解法及其各种变形。 4.教学方法:讲授与讨论。 二、教学过程 在上节中我们用矩阵初等变换来分析Gauss 消去法,得到了重要的矩阵LU 分解定理(定理 3.1,3.2)。由此我们将得到Gauss 消去法的变形:直接三角分解法。直接三角分解法的基本想法是,一旦实现了矩阵A 的LU 分解,那么求解方程组b x =A 的问题就等价于求解两个三角形方程组 (1)b y =L ,求y ; (2)y x =U ,求x 。 而这两个三角形方程组的求解是容易的。下面我们先给出这两个三角形方程组的求解公式;然后研究在LU A =或LU PA =时,U L ,的元素与A 的元素之间的直接关系。 4-0 三角形线性方程组的解法 设 ????? ???????= nn n n l l l l l l L 21222111, 11121222n n nn u u u u u U u ??????=???????? 则b y =L 为下三角形方程组,它的第i 个方程为 ),2,1(11,22111 n i b y l y l y l y l y l i i ii i i i i i i j j ij ==++++=--=∑ 假定0≠ii l ,按n y y y ,,,21 的顺序解得: ??? ?? ? ?=+-==∑-=) ,,3,2(/1111 11n i l b y l y l b y ii i i j j ij i 上三角形方程组y x =U 的第i 个方程为

解三角形应用举例练习高考试题练习

解三角形应用举例练习 班级 姓名 学号 得分 一、选择题 1.从A 处望B 处的仰角为α,从B 处望A 处的俯角为β,则α、β的关系为…………………( ) A.α>β B.α=β C.α+β=90° D.α+β=180° 2.在200米高的山顶上,测得山下一塔顶与塔底的俯角分别为30°、60°,则塔高为…..( ) A. 3 400 B. 33400米 C. 2003米 D. 200米 3.在?ABC 中, 已知sinA = 2 sinBcosC, 则?ABC 一定是…………………………………….( ) A. 直角三角形; B. 等腰三角形; C.等边三角形; D.等腰直角三角形. 4.如图,△ABC 是简易遮阳棚,A 、B 是南北方向上两个定点,正东方向射出的太阳光线与地面 成40°角,为了使遮阴影面ABD 面积最大,遮阳棚ABC 与地面所成的角为……………….( ) A C D B 阳光地面 A.75° B.60° C.50° D.45° 5.台风中心从A 地以20 km/h 的速度向东北方向移动,离台风中心30 km 内的地区为危险区,城市B 在A 的正东40 km 处,B 城市处于危险区内的时间为…………………………………..( ) A.0.5 h B.1 h C.1.5 h D.2 h 6.在△ABC 中,已知b = 6,c = 10,B = 30°,则解此三角形的结果是 …………………( ) A 、无解 B 、一解 C 、两解 D 、解的个数不能确定 二、填空题 7. 甲、乙两楼相距20米,从乙楼底望甲楼顶的仰角为60°,从甲楼顶望乙楼顶的俯角为30°,则甲、乙两楼的高分别是 8.我舰在敌岛A 南50°西相距12nmile 的B 处,发现敌舰正由岛沿北10°西的方向以10nmile/h 的速度航行,我舰要用2小时追上敌舰,则需要速度的大小为 9.有一两岸平行的河流,水速为1,小船的速度为2,为使所走路程最短,小船应朝_______方 向行驶. C D 12 A B D 6045 0 m o o 10..在一座20 m 高的观测台顶测得地面一水塔塔顶仰角为60°,塔底俯角为45°,那么这座塔的 高为_______.

列主元三角分解法在matlab中的实现

列主元三角分解法在matlab中的实现 摘要:介绍了M atlab语言并给出用M atlab语言实现线性方程组的列主元三角分解法,其有效性已在计算机实现中得到了验证。 关键词:M atlab语言;高斯消去法;列主元三角分解法 0前言 M atlab是M atrix Laboratory(矩阵实验室)的缩写,它是由美国M athwork公司于1967年推出的软件包,现已发展成为一种功能强大的计算机语言。它编程简单,使用方便,在M a tlab环境下数组的操作与数的操作一样简单,进行数学运算可以像草稿纸一样随心所欲,使计算机兼备高级计算器的优点。M atlab语言具有强大的矩阵和向量的操作功能,是Fo rtran和C语言无法比拟的;M a tlab语言的函数库可任意扩充;语句简单,内涵丰富;还具有二维和三维绘图功能且使用方便,特别适用于科学和工程计算。 在科学和工程计算中,应用最广泛的是求解线性方程组的解,一般可用高斯消去法求解,如果系数矩阵不满足高斯消去法在计算机上可行的条件,那么消元过程中可能会出现零主元或小主元,消元或不可行或数值不稳定,解决办法就是对方程组进行行交换或列交换来消除零主元或小主元,这就是选主元的思想。 1 定义 列主元三角分解:如果A为非奇异矩阵,则存在排列矩阵P,使PA=LU,其中L为单位下三角矩阵,U为上三角阵。列主元三角分角法是对直接三角分解法的一种改进,主要目的和列主元高斯消元法一样,

就是避免小数作为分母项. 2 算法概述 列主元三角分解法和普通三角分解法基本上类似,所不同的是在构造Gauss 变换前,先在对应列中选择绝对值最大的元素(称为列主元),然后实施初等行交换将该元素调整到矩阵对角线上。 例如第)1,,2,1(-=n k 步变换叙述如下: 选主元:确定p 使{}1)1( max -≤≤-=k ik n i k k pk a a ; 行交换:将矩阵的第k 行和第p 行上的元素互换位置,即 . 实施Gauss 变换:通过初行变换,将列主对角线以下的元素消为零.即 3 列主元三角分解在matlab 中的实现

最新解三角形应用举例练习题

解三角形应用举例练习题 一、选择题 1.某人向正东方向走x km后,他向右转150°,然后朝新方向走3 km,结果他离出发点恰好 3 km,那么x的值为() A.3B.2 3 C.23或 3 D.3 2.已知船A在灯塔C北偏东85°且到C的距离为2km,船B在灯塔C西偏北25°且到C的距离为3km,则A,B两船的距离为() A.23km B.32km C.15km D.13km 3.已知△ABC的三边长a=3,b=5,c=6,则△ABC的面积是() A.14 B.214 C.15 D.215 4.两座灯塔A和B与海洋观察站C的距离都等于a km,灯塔A在观察站C的北偏东20°,灯塔B在观察站C的南偏东40°,则灯塔A与灯塔B的距离为() A.a km B.3a km C.2a km D.2a km 5.已知△ABC中,a=2、b=3、B=60°,那么角A等于() A.135°B.90° C.45°D.30° 6.一船向正北航行,看见正西方向有相距10海里的两个灯塔恰好与它在一条直线上,继续航行半小时后,看见一灯塔在船的南偏西60°方向上,另一灯塔在船的南偏西75°方向上,则这艘船的速度是每小时() A.5海里B.53海里 C.10海里D.103海里 二、填空题 7.(2010~2011·醴陵二中、四中期中)已知A、B两地的距离为10km,BC两地的距离

为20km,经测量∠ABC=120°,则AC两地的距离为________km. 8.如图,为了测量河的宽度,在一岸边选定两点A,B,望对岸的标记物C,测得∠CAB=30°,∠CBA=75°,AB=120 m,则河的宽度是__________. 9. (2011·北京朝阳二模)如图,一艘船上午在A处测得灯塔S在它的北偏东30°处,之后它继续沿正北方向匀速航行,上午到达B处,此时又测得灯塔S在它的北偏东75°处,且与它相距42n mile,则此船的航行速度是________n mile/h. 三、解答题

(完整版)有限元大作业matlab---课程设计例子

有限元大作业程序设计 学校:天津大学 院系:建筑工程与力学学院 专业:01级工程力学 姓名:刘秀 学号:\\\\\\\\\\\ 指导老师:

连续体平面问题的有限元程序分析 [题目]: 如图所示的正方形薄板四周受均匀载荷的作用,该结构在边界 上受正向分布压力, m kN p 1=,同时在沿对角线y 轴上受一对集中压 力,载荷为2KN ,若取板厚1=t ,泊松比0=v 。 [分析过程]: 由于连续平板的对称性,只需要取其在第一象限的四分之一部分参加分析,然后人为作出一些辅助线将平板“分割”成若干部分,再为每个部分选择分析单元。采用将此模型化分为4个全等的直角三角型单元。利用其对称性,四分之一部分的边界约束,载荷可等效如图所示。

[程序原理及实现]: 用FORTRAN程序的实现。由节点信息文件NODE.IN和单元信息文件ELEMENT.IN,经过计算分析后输出一个一般性的文件DATA.OUT。模型基本信息由文件为BASIC.IN生成。 该程序的特点如下: 问题类型:可用于计算弹性力学平面问题和平面应变问题 单元类型:采用常应变三角形单元 位移模式:用用线性位移模式 载荷类型:节点载荷,非节点载荷应先换算为等效节点载荷 材料性质:弹性体由单一的均匀材料组成 约束方式:为“0”位移固定约束,为保证无刚体位移,弹性体至少应有对三个自由度的独立约束 方程求解:针对半带宽刚度方程的Gauss消元法

输入文件:由手工生成节点信息文件NODE.IN,和单元信息文件ELEMENT.IN 结果文件:输出一般的结果文件DATA.OUT 程序的原理如框图:

选主元的三角分解法

2012-2013(1)专业课程实践论文选主元(列)的三角分解法 范俊,0818180124,R数学08-1班 毛勇,0818180117,R数学08-1班

一、算法理论 从直接三角分解公式可看出,当时计算将中断或者当绝对值很小时,按分解公式计算可能要求舍入误差的累积,但如果非奇异,就可通过交换的行实现矩阵的分解,因此可采用与列主元素消去法类似的方法 (可以证明下述方法与列主元素消去法等价),将直接三角分解法修改为(部分)选主元的三角分解法。 设第步分解已完成,这时有 第步分解需要用到式 及式 为了避免用小的数作除数,引进量 于是有 ,,

用作为,交换的行与行元素(将位置的新元素仍记作及 ),于是有。由此在进行第步分解计算。 该程序即利用选主元的三角分解法解方程求方程的根。先选择列主元,再构造LU矩阵,通过求解LY = B和UX = Y得出需要的解。 注:方程维数在程序中需按题意自定义。

源代码源代码: LU_Decomposition.cpp #include #include #define N 4 //矩阵维数,可自定义 static double A[N][N]; //系数矩阵 static double B[N]; //右端项 static double Y[N]; //中间项 static double X[N]; //输出 static double S[N]; //选取列主元的比较器 int i,j,k; //计数器 void main() { cout << "请输入线性方程组(ai1,ai2,ai3......ain, yi):"<> A[i][j]; cin >>B[i]; } for (k = 0 ; k < N; k++) { //选列主元 int index =k; for (i = k ; i < N; i++) { double temp = 0; for (int m = 0 ; m < k ;m++) { temp = temp + A[i][m]*A[m][k]; } S[i] = A[i][k] - temp; if(S[index] < S[i]) { index = i; }

解三角形应用举例

东方中学教案 1.知识与技能: 会在各种应用问题中,抽象或构造出三角形,标出已知量、未知量,确定解三角形的方法;搞清利用解斜三角形可解决的各类应用问题的基本图形和基本等量关系;理解各种应用问题中的有关名词、术语,如:坡度、俯角、仰角、方向角、方位角等;通过解三角形的应用的学习,提高解决实际问题的能力 2.过程与方法: 通过巧妙的设疑,顺利的引导新课,为下节课做好铺垫。结合学生的实际情况,采用“提出问题—引发思考—探索猜想—总结规律—反馈练习”的教学过程,根据大纲要求以及教学内容之间的内在联系,铺开例题,设计变式,同时通过多媒体、图形观察等直观演示,帮助学生掌握在各种应用问题中,抽象或构造出三角形,标出已知量、未知量,确定解三角形的方法。 3.情感、态度与价值观: 实际问题中抽象出一个或几个三角形,然后逐个解三角形,得到实际问题的解。

修改简记教学过程: 一、复习引入: 二、讲解范例: 例1 自动卸货汽车的车箱采用液压结构,设计时需要计算 油泵顶杆BC的长度已知车箱的最大仰角为60°,油泵顶点 B与车箱支点A之间的距离为1.95m,AB与水平线之间的夹角 为6°20′,AC长为1.40m,计算BC的长(保留三个有效数字) 分析:求油泵顶杆BC的长度也就是在△ABC内,求边长BC的问题,而根据已知条件, AC=1.40m,AB=1.95 m,∠BAC=60°+6°20′=66°20′相当于已知△ABC 的两边和它们的夹角,所以求解BC可根据余弦定理解:由余弦定理,得 BC2=AB2+AC2-2AB·AC cos A =1.952+1.402-2×1.95×1.40×cos66°20′=3.571 ∴BC≈1.89 (m) 答:油泵顶杆B C约长1.89 m 评述:此题虽为解三角形问题的简单应用,但关键是把未知边所处的三角形找到,在转 换过程中应注意“仰角”这一概念的意义,并排除题目中非数学因素的干扰,将数量关系 从题目准确地提炼出来 例2某渔船在航行中不幸遇险,发出求救信号,我海军舰艇在A处获悉后,立即测出该渔 船在方位角为45°、距离A为10海里的C处,并测得渔船正沿方位角为105°的方向, 以9海里/h的速度向某小岛B靠拢,我海军舰艇立即以21海里/h的速度前去营救, 试问舰艇应按照怎样的航向前进?并求出靠近渔船所用的时间

Matlab有限元分析操作基础共11页

Matlab有限元分析20140226 为了用Matlab进行有限元分析,首先要学会Matlab基本操作,还要学会使用Matlab进行有限元分析的基本操作。 1. 复习:上节课分析了弹簧系统 x 推导了系统刚度矩阵

2. Matlab有限元分析的基本操作 (1)单元划分(选择何种单元,分成多少个单元,标号)(2)构造单元刚度矩阵(列出…) (3)组装系统刚度矩阵(集成整体刚度矩阵) (4)引入边界条件(消除冗余方程) (5)解方程 (6)后处理(扩展计算)

3. Matlab有限元分析实战【实例1】

分析: 步骤一:单元划分

>>k1=SpringElementStiffness(100)

a) 分析SpringAssemble库函数 function y = SpringAssemble(K,k,i,j) % This function assembles the element stiffness % matrix k of the spring with nodes i and j into the % global stiffness matrix K. % function returns the global stiffness matrix K % after the element stiffness matrix k is assembled. K(i,i) = K(i,i) + k(1,1); K(i,j) = K(i,j) + k(1,2); K(j,i) = K(j,i) + k(2,1); K(j,j) = K(j,j) + k(2,2); y = K; b) K是多大矩阵? 今天的系统刚度矩阵是什么? 因为 11 22 1212 k k k k k k k k - ?? ?? - ????--+ ?? 所以 1000100 0200200 100200300 - ?? ?? - ????-- ???

解三角形应用举例最新衡水中学自用精品教学设计

解三角形应用举例 主标题:解三角形应用举例 副标题:为学生详细的分析解三角形应用举例的高考考点、命题方向以及规律总结。 关键词:距离测量,高度测量,仰角,俯角,方位角,方向角 难度:3 重要程度:5 考点剖析: 能够运用正弦定理、余弦定理等知识解决一些与测量和几何计算有关的实际问题. 命题方向: 1.测量距离问题是高考的常考内容,既有选择、填空题,也有解答题,难度适中,属中档题. 2.高考对此类问题的考查常有以下两个命题角度: (1)测量问题; (2)行程问题. 规律总结: 1个步骤——解三角形应用题的一般步骤 2种情形——解三角形应用题的两种情形 (1)实际问题经抽象概括后,已知量与未知量全部集中在一个三角形中,可用正弦定理或余弦定理求解. (2)实际问题经抽象概括后,已知量与未知量涉及到两个或两个以上的三角形,这时需作出这些三角形,先解够条件的三角形,然后逐步求解其他三角形,有时需设出未知量,从几个三角形中列出方程(组),解方程(组)得出所要求的解. 2个注意点——解三角形应用题应注意的问题 (1)画出示意图后要注意寻找一些特殊三角形,如等边三角形、直角三角形、等腰三角形等,这样可以优化解题过程. (2)解三角形时,为避免误差的积累,应尽可能用已知的数据(原始数据),少用间接求出的量.

知识梳理 1.距离的测量 背景可测元素图形目标及解法 两点均可到达a,b,α 求AB:AB= a2+b2-2ab cos α 只有一点可到达b,α,β 求AB:(1)α+β+B=π; (2) AB sin β= b sin B 两点都不可到达a,α,β, γ,θ 求AB:(1)△ACD中,用 正弦定理求AC; (2)△BCD中,用正弦定理 求BC; (3)△ABC中,用余弦定理 求AB 2.高度的测量 背景可测元素图形目标及解法 底部可 到达 a,α求AB:AB=a tan_α 底部不可到达a,α,β 求AB:(1)在△ACD中用正弦 定理求AD;(2)AB=AD sin_β 3.实际问题中常见的角 (1)仰角和俯角 在同一铅垂平面内的水平视线和目标视线的夹角,目标视线在水平视线上方时叫仰角,目标视线在水平视线下方时叫俯角(如图1).

基于matlab的有限元法分析平面应力应变问题刘刚

姓名:刘刚学号:15 平面应力应变分析有限元法 Abstruct:本文通过对平面应力/应变问题的简要理论阐述,使读者对要分析的问题有大致的印象,然后结合两个实例,通过MATLAB软件的计算,将有限元分析平面应力/应变问题的过程形象的展示给读者,让人一目了然,快速了解有限元解决这类问题的方法和步骤! 一.基本理论 有限元法的基本思路和基本原则以结构力学中的位移法为基础,把复杂的结构或连续体看成有限个单元的组合,各单元彼此在节点出连接而组成整体。把连续体分成有限个单元和节点,称为离散化。先对单元进行特性分析,然后根据节点处的平衡和协调条件建立方程,综合后做整体分析。这样一分一合,先离散再综合的过程,就是把复杂结构或连续体的计算问题转化简单单元分析与综合问题。因此,一般的有限揭发包括三个主要步骤:离散化单元分析整体分析。 二.用到的函数 1. LinearTriangleElementStiffness(E,NU,t,xi,yi,xj,yj,xm,ym,p) (K k I f) (k u) (k u A) (E NU t) 三.实例 例1.考虑如图所示的受均布载荷作用的薄平板结构。将平板离散化成两个线性三角元,假定E=200GPa,v=,t=0.025m,w=3000kN/m. 1.离散化 2.写出单元刚度矩阵

通过matlab 的LinearTriangleElementStiffness 函数,得到两个单元刚度矩阵1k 和2k ,每个矩阵都是6 6的。 >> E=210e6 E = >> k1=LinearTriangleElementStiffness(E,NU,t,0,0,,,0,,1) k1 = +006 * Columns 1 through 5 0 0 0 0 0 0 0 0 Column 6 >> NU= NU = >> t= t = >> k2=LinearTriangleElementStiffness(E,NU,t,0,0,,0,,,1)

列主元高斯消去法和列主元三角分解法解线性方程

列主元高斯消去法和列主元三角分解法解线性方程

计算方法实验报告1 【课题名称】 用列主元高斯消去法和列主元三角分解法解线性方程 【目的和意义】 高斯消去法是一个古老的求解线性方程组的方法,但由它改进得到的选主元的高斯消去法则是目前计算机上常用的解低阶稠密矩阵方程组的有效方法。 用高斯消去法解线性方程组的基本思想时用矩阵行的初等变换将系数矩阵A 约化为具有简单形式的矩阵(上三角矩阵、单位矩阵等),而三角形方程组则可以直接回带求解 用高斯消去法解线性方程组b Ax =(其中A ∈Rn ×n )的计算量为:乘除法运算步骤为 32(1)(1)(21)(1)(1)262233 n n n n n n n n n n n MD n ----+= +++=+-,加减运算步骤为 (1)(21)(1)(1)(1)(25) 6226 n n n n n n n n n n AS -----+= ++= 。相比之 下,传统的克莱姆法则则较为繁琐,如求解20阶线性方程组,克莱姆法则大约要19 510?次乘法, 而用高斯消去法只需要3060次乘除法。

在高斯消去法运算的过程中,如果出现abs(A(i,i))等于零或过小的情况,则会导致矩阵元素数量级严重增长和舍入误差的扩散,使得最后的计算结果不可靠,所以目前计算机上常用的解低阶稠密矩阵方程的快速有效的方法时列主元高斯消去法,从而使计算结果更加精确。 2、列主元三角分解法 高斯消去法的消去过程,实质上是将A 分解为两个三角矩阵的乘积A=LU ,并求解Ly=b 的过程。回带过程就是求解上三角方程组Ux=y 。所以在实际的运算中,矩阵L 和U 可以直接计算出,而不需要任何中间步骤,从而在计算过程中将高斯消去法的步骤进行了进一步的简略,大大提高了运算速度,这就是三角分解法 采用选主元的方式与列主元高斯消去法一样,也是为了避免除数过小,从而保证了计算的精确度 【计算公式】 1、 列主元高斯消去法 设有线性方程组Ax=b ,其中设A 为非奇异矩阵。方程组的增广矩阵为 1112112122221[,]n n l a a a a b a a a b ?????? ??=? ? ?? ??? ?a b L L M L L M M M M M M M M M

解三角形应用举例

第三章 三角函数、三角恒等变换及解三角形第8课时 解三角 形应用举例 1. (必修5P 11习题4改编)若海上有A 、B 、C 三个小岛,测得A ,B 两岛相距10海里,∠BAC =60°,∠ABC =75°,则B 、C 间的距离是________海里. 答案:5 6 解析:由正弦定理, 知 BC sin60°=AB sin (180°-60°-75°) , 解得BC =56(海里). 2. (必修5P 20练习第4题改编)江岸边有一炮台高30 m ,江中有两条船,船与炮台底部在同一水面上,由炮台顶部测得俯角分别为45°和60°,而且两条船与炮台底部连线成30°角,则两条船相距________m. 答案:10 3 解析:如图,OA 为炮台,M 、N 为两条船的位置,∠AMO =45°,∠ANO =60°,OM =AOtan45°=30,ON =AOtan30°= 3 3 ×30=103,由余弦定理,得 MN = 900+300-2×30×103× 3 2 =300=103(m). 3. (必修5P 18例1改编)如图,要测量河对岸A 、B 两点间的距离,今沿河岸选取相距40 m 的C 、D 两点,测得∠ACB=60°,∠BCD =45°,∠ADB =60°,∠ADC =30°,则AB 的距离是__________ m. 答案:20 6 解析:由已知知△BDC 为等腰直角三角形,故DB =40;由∠ACB=60°和∠ADB=60°知A 、B 、C 、D 四点共圆, 所以∠BAD=∠BCD=45°;

在△BDA 中,运用正弦定理可得AB =20 6. 4. (必修5P 21习题2改编)某人在C 点测得塔顶A 在南偏西80°,仰角为45°,此人沿南偏东40°方向前进10 m 到D ,测得塔顶A 的仰角为30°,则塔高为________m. 答案:10 解析:如图,设塔高为h ,在Rt △AOC 中,∠ACO =45°,则OC =OA =h. 在Rt △AOD 中,∠ADO =30°,则OD =3h. 在△OCD 中,∠OCD =120°,CD =10. 由余弦定理得OD 2=OC 2+CD 2 -2OC·CD cos ∠OCD , 即(3h)2 =h 2 +102 -2h×10×cos120°, ∴ h 2 -5h -50=0,解得h =10或h =-5(舍). 5. 如图,一船在海上自西向东航行,在A 处测得某岛M 的方位角为北偏东α角,前进mkm 后在B 处测得该岛的方位角为北偏东β角,已知该岛周围nkm 范围内(包括边界)有暗礁,现该船继续东行.当α与β满足条件________时,该船没有触礁危险. 答案:mcos αcos β>nsin(α-β) 解析:∠MAB=90°-α,∠MBC =90°-β=∠MAB+∠AMB=90°-α+∠AMB,∴ ∠AMB =α-β.由题可知,在△ABM 中,根据正弦定理得BM sin (90°-α)=m sin (α-β), 解得BM = mcos αsin (α-β).要使船没有触礁危险,需要BMsin(90°-β)=mcos αcos β sin (α-β) >n , 所以α与β满足mcos αcos β>nsin(α-β)时船没有触礁危险. 1. 用正弦定理和余弦定理解三角形的常见题型 测量距离问题、高度问题、角度问题、计算面积问题、航海问题、物理问题等. 2. 实际问题中的常用角 (1) 仰角和俯角 与目标线在同一铅垂平面内的水平视线和目标视线的夹角,目标视线在水平视线上方的角叫仰角,目标视线在水平视线下方的角叫俯角(如图①). (2) 方向角:相对于某正方向的水平角,如南偏东30°,北偏西45°,西偏北60°等.

解三角形应用举例

第7节 解三角形应用举例 最新考纲 能够运用正弦定理、余弦定理等知识方法解决一些与测量、几何计算有关的实际问题. 知 识 梳 理 1.仰角和俯角 在同一铅垂平面内的水平视线和目标视线的夹角,目标视线在水平视线上方叫仰角,目标视线在水平视线下方叫俯角(如图1). 2.方向角 相对于某正方向的水平角,如南偏东30°,北偏西45°等. 3.方位角 指从正北方向顺时针转到目标方向线的水平角,如B 点的方位角为α(如图2). 4.坡度:坡面与水平面所成的二面角的正切值. [常用结论与微点提醒] 1.不要搞错各种角的含义,不要把这些角和三角形内角之间的关系弄混. 2.在实际问题中,可能会遇到空间与平面(地面)同时研究的问题,这时最好画两个图形,一个空间图形,一个平面图形,这样处理起来既清楚又不容易出现错误. 诊 断 自 测 1.思考辨析(在括号内打“√”或“×”) (1)东北方向就是北偏东45°的方向.( ) (2)从A 处望B 处的仰角为α,从B 处望A 处的俯角为β,则α,β的关系为α+β=180°.( ) (3)俯角是铅垂线与视线所成的角,其范围为? ?????0,π2.( ) (4)方位角与方向角其实质是一样的,均是确定观察点与目标点之间的位置关系.( )

解析 (2)α=β;(3)俯角是视线与水平线所构成的角. 答案 (1)√ (2)× (3)× (4)√ 2.若点A 在点C 的北偏东30°,点B 在点C 的南偏东60°,且AC =BC ,则点A 在点B 的( ) A.北偏东15° B.北偏西15° C.北偏东10° D.北偏西10° 解析 如图所示,∠ACB =90°, 又AC =BC , ∴∠CBA =45°,而β=30°, ∴α=90°-45°-30°=15°. ∴点A 在点B 的北偏西15°. 答案 B 3.(教材习题改编)如图所示,设A ,B 两点在河的两岸,一测量 者在A 所在的同侧河岸边选定一点C ,测出AC 的距离为50 m , ∠ACB =45°,∠CAB =105°后,就可以计算出A ,B 两点的 距离为( ) A.50 2 m B.50 3 m C.25 2 m D.2522 m 解析 由正弦定理得AB sin ∠ACB =AC sin B , 又∵B =30°,∴AB =AC sin ∠ACB sin B =50×2212 =502(m). 答案 A 4.轮船A 和轮船B 在中午12时同时离开海港C ,两船航行方向的夹角为120°,两船的航行速度分别为25 n mile/h ,15 n mile/h ,则下午2时两船之间的距离是______n mile. 解析 设两船之间的距离为d , 则d 2=502+302-2×50×30×cos 120°=4 900, ∴d =70,即两船相距70 n mile.

PDE的Galerkin和有限元的MATLAB程序

%Galerkin方法 clear all clc syms x; A=zeros(4,4); for i=1:4 for j=1:4 phy1=x^i; phy2=x^j; dphy1=diff(phy1,1); dphy2=diff(phy2,1); phy=pi^2*phy1*phy2; dphy=dphy1*dphy2; A(i,j)=int(phy+dphy,x,0,1); end end D=[]; for k=1:4 f1=2*pi^2*sin(pi*x)+pi^3*x; f2=x^k; f=f1*f2; D(k)=int(f,x,0,1); end D=D'; C=A\D; C=C'; X=linspace(0,1,6); F=0; for i=1:4 F=F+C(i)*x^i; end for j=1:6 Y(j)=subs(F,X(j)); end Y=Y-pi.*X; Y1=sin(pi.*X); err=norm(abs(Y-Y1)); disp('数值解') disp(Y) disp('整体误差') disp(err)

%%%%%%%%%追赶法 function x=chase(a,b,c,d) n=length(b); u(1)=c(1)/b(1); q(1)=d(1)/b(1); for i=2:n-1 h(i)=b(i)-u(i-1)*a(i-1); u(i)=c(i)/h(i); q(i)=(d(i)-q(i-1)*a(i-1))/h(i); end q(n)=(d(n)-q(n-1)*a(n-1))/(b(n)-u(n-1)*a(n-1)); x(n)=q(n); for i=n-1:-1:1 x(i)=q(i)-u(i)*x(i+1); end end %有限元法 clear all clc m=5; xspan=[0 1]; a=xspan(1);b=xspan(2);h=(b-a)/m; x=[a:h:b]; f=zeros(m-1,1); F1=@(t)(-1/h+(pi^2)*t.*(1-t)*h); c1=quadl(F1,0,1); F2=@(t)(2/h+(pi^2)*(t.^2+(1-t).^2)*h); c2=quadl(F2,0,1); A=diag(c2*ones(m-1,1))+diag(c1*ones(m-2,1),1)+... diag(c1*ones(m-2,1),-1); for i=2:m F=@(t)h*t.*(2*pi^2*sin(pi*(x(i-1)+h*t)))+... h*(1-t).*(2*pi^2*sin(pi*(x(i)+h*t))); f(i-1)=quadl(F,0,1); end v=chase(c1*ones(m-2,1),c2*ones(m-1,1),c1*ones(m-2,1),f); u=[0,v,0]

matlab有限元分析实例

1.物理现象:这个对工程师来说是直观的物理现象和物理量,温 度多少度,载荷是多大等等。通常来说,用户界面中呈现的、用户对工程问题进行设置时输入的都是此类信息。 2.数学方程:将物理现象翻译成相应的数学方程,例如流体对应 的是NS方程,传热对应的是传热方程等等;大部分描述这些现象的方程在空间上都是偏微分方程,偶尔也有ODE(如粒子轨迹、化学反应等)。在这个层面,软件把物理现象“翻译” 为以解析式表示的数学模型。 3.数值模型:在定义了数学模型,并执行了网格剖分后,商业软 件会将数学模型离散化,利用有限元方法、边界元法、有限差分法、不连续伽辽金法等方法生成数值模型。软件会组装并计算方程组雅可比矩阵,并利用求解器求解方程组。这个层面的计算通常是隐藏在后台的,用户只能通过一些求解器的参数来干预求解。 有限元是一种数值求解偏微分方程的方法。 基本过程大致是设置形函数,离散,形成求解矩阵,数值解矩阵,后处理之类的。 MATLAB要把这些过程均自己实现,不过在数值求解矩阵时可以调用已有函数。可以理解为MATLAB是一个通用的计算器,当然它的功能远不止如此。

而ANSYS之类的叫做通用有限元软件,针对不同行业已经将上述过程封装,前后处理也比较漂亮,甚至不太了解有限元理论的人也能算些简单的东西,当然结果可靠性又另说了。 比较两者,ANSYS之类的用起来容易得多,但灵活性不如MATLAB。MATLAB用起来很困难,也有人做了一些模块,但大多数只能解决一些相对简单的问题。 对于大多数工程问题,以及某些领域的物理问题,一般都用通用有限元软件,这些软件还能添加一些函数块,用以解决一些需要额外设置的东西。但是对于非常特殊的问题,以及一般性方程的有限元解,那只能用MATLAB或C,Fortran之类的了。

有限元的MATLAB解法

有限元的MATLAB解法 1.打开MATLAB。 2.输入“pdetool”再回车,会跳出PDE Toolbox的窗口(PDE意为偏微分方程,是partial differential equations的缩写),需要的话可点击Options菜单下Grid命令,打开栅格。 3.完成平面几何模型:在PDE Toolbox的窗口中,点击工具栏下的矩形几何模型进行制作模型,可画矩形R,椭圆E,圆C,然后在Set formula栏进行编辑并(如双脊波导R1+R2+R3改为RI-R2-R3,设定a、b、s/a、d/b的值从而方便下步设定坐标) 用算术运算符将图形对象名称连接起来,若还需要,可进行储存,形成M文件。 4.用左键双击矩形进行坐标设置:将大的矩形left和bottom都设为0,width是矩形波导的X轴的长度,height是矩形波导的y轴的长度,以大的矩形左下角点为原点坐标为参考设置其他矩形坐标。 5.进行边界设置:点击“Boundary”中的“Boundary Mode”,再点

击“Boundary”中的“Specify Boundary Conditions”,选择符合的边界条件,Neumann为诺曼条件,Dirichlet为狄利克雷条件,边界颜色显示为红色。 6.进入PDE模式:点击"PDE"菜单下“PDE Mode”命令,进入PDE模式,单击“PDE Specification”,设置方程类型,“Elliptic”为椭圆型,“Parabolic”为抛物型,“Hyperbolic”为双曲型,“Eigenmodes”为特征值问题。 7.对模型进行剖分:点击“Mesh”中“Initialize Mesh”进行初次剖分,若要剖的更细,再点击“Refine Mesh”进行网格加密。 8.进行计算:点击“Solve”中“Solve PDE”,解偏微分方程并显示图形解,u值即为Hz或者Ez。 9.单击“Plot”菜单下“Parameters”选项,打开“Plot Selection”对话框。选中Color,Height(3-D plot)和Show mesh三项,然后单击“Plot”按钮,显示三维图形解。 10.如果要画等值线图和矢量场图,单击“Plot”菜单下“Parameters”选项,打开“Plot Selection”对话框。选中Contour和Arrows两项,然后单击Plot按钮,可显示解的等值线图和矢量场图。 11.将计算结果条件和边界导入MATLAB中:点击“Export Solution”,再点击“Mesh”中“Export Mesh”。

相关主题
文本预览
相关文档 最新文档