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有限元分析与应用》由简单到复杂，循序渐进地介绍了各种有限元及其分析与应用方法。书中提供了大量取自机械工程、土木工程、航空航天工程和材料科学的示例和习题，具有很高的工程应用价值。
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有限元分析与应用》由简单到复杂，循序渐进地介绍了各种有限元及其分析与应用方法。书中提供了大量取自机械工程、土木工程、航空航天工程和材料科学的示例和习题，具有很高的工程应用价值。
数值分析MATLAB上机实验
数值分析实习报告 姓名：gestepoA 学号：201******* 班级：***班
序言 随着计算机技术的迅速发展，数值分析在工程技术领域中的应用越来越广泛，并且成为数学与计算机之间的桥梁。要解决工程问题，往往需要处理很多数学模型，不仅要研究各种数学问题的数值解法，同时也要分析所用的数值解法在理论上的合理性，如解法所产生的误差能否满足精度要求：解法是否稳定、是否收敛及熟练的速度等。而且还能减少大量的人工计算。 由于工程实际中所遇到的数学模型求解过程迭代次数很多，计算量很大，所以需要借助如MATLAB，C++，VB，JAVA的辅助软件来解决，得到一个满足误差限的解。本文所计算题目，均采用MATLAB进行编程，MATLAB被称为第四代计算机语言，利用其丰富的函数资源，使编程人员从繁琐的程序代码中解放出来MATLAB最突出的特点就是简洁，它用更直观的、符合人们思维习惯的代码。它具有以下优点： 1友好的工作平台和编程环境。MATLAB界面精致，人机交互性强，操作简单。 2简单易用的程序语言。MATLAB是一个高级的矩阵/阵列语言，包含控制语言、函数、数据结构，具有输入、输出和面向对象编程特点。用户可以在命令窗口中将输入语句与执行命令同步，也可以先编好一个较大的复杂的应用程序（M 文件）后再一起运行。 3强大的科学计算机数据处理能力。包含大量计算算法的集合，拥有600多个工程中要用到的数学运算函数。 4出色的图像处理功能，可以方便地输出二维图像，便于我们绘制函数图像。
目录 1 第一题 (4) 1.1 实验目的 (4) 1.2 实验原理和方法 (4) 1.3 实验结果 (5) 1.3.1 最佳平方逼近法 (5) 1.3.2 拉格朗日插值法 (7) 1.3.3 对比 (8) 2 第二题 (9) 2.1实验目的 (9) 2.2 实验原理和方法 (10) 2.3 实验结果 (10) 2.3.1 第一问 (10) 2.3.2 第二问 (11) 2.3.3 第三问 (11) 3 第三题 (12) 3.1实验目的 (12) 3.2 实验原理和方法 (12) 3.3 实验结果 (12) 4 MATLAB程序 (14)
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 - ?? ?? - ?? ?? -- ?? ？
《MATLAB与数值分析》第一次上机实验报告
电子科技大学电子工程学院标准实验报告（实验）课程名称MATLAB与数值分析 学生姓名：李培睿 学号：2013020904026 指导教师：程建
一、实验名称 《MATLAB与数值分析》第一次上机实验 二、实验目的 1. 熟练掌握矩阵的生成、加、减、乘、除、转置、行列式、逆、范数等运算 操作。（用.m文件和Matlab函数编写一个对给定矩阵进行运算操作的程序） 2. 熟练掌握算术符号操作和基本运算操作，包括矩阵合并、向量合并、符号 转换、展开符号表达式、符号因式分解、符号表达式的化简、代数方程的符号解析解、特征多项式、函数的反函数、函数计算器、微积分、常微分方程的符号解、符号函数的画图等。（用.m文件编写进行符号因式分解和函数求反的程序） 3. 掌握Matlab函数的编写规范。 4、掌握Matlab常用的绘图处理操作，包括：基本平面图、图形注释命令、 三维曲线和面的填充、三维等高线等。（用.m文件编写在一个图形窗口上绘制正弦和余弦函数的图形，并给出充分的图形注释） 5. 熟练操作MATLAB软件平台，能利用M文件完成MATLAB的程序设计。 三、实验内容 1. 编程实现以下数列的图像，用户能输入不同的初始值以及系数。并以x， y为坐标显示图像 x(n+1) = a*x(n)-b*(y(n)-x(n)^2); y(n+1) = b*x(n)+a*(y(n)-x(n)^2) 2. 编程实现奥运5环图，允许用户输入环的直径。 3. 实现对输入任意长度向量元素的冒泡排序的升序排列。不允许使用sort 函数。 四、实验数据及结果分析 题目一： ①在Editor窗口编写函数代码如下：
数值分析的matlab实现
第2章牛顿插值法实现 参考文献：[1]岑宝俊. 牛顿插值法在凸轮曲线修正设计中的应用[J]. 机械工程师,2009,10:54-55. 求牛顿插值多项式和差商的MA TLAB 主程序： function[A,C,L,wcgs,Cw]=newpoly(X,Y) n=length(X);A=zeros(n,n);A(:,1) =Y'; s=0.0;p=1.0;q=1.0;c1=1.0; for j=2:n for i=j:n A(i,j)=(A(i,j-1)-A(i-1,j-1))/(X(i)-X(i-j+1)); end b=poly(X(j-1));q1=conv(q,b);c1=c1*j;q=q1; end C=A(n,n);b=poly(X(n));q1=conv(q1,b); for k=(n-1):-1:1 C=conv(C,poly(X(k)));d=length(C);C(d)=C(d)+A(k,k); end L(k,:)=poly2sym(C);Q=poly2sym(q1); syms M wcgs=M*Q/c1;Cw=q1/c1; （1）保存名为newpoly.m 的M 文件 （2）输入MA TLAB 程序 >> X=[242,243,249,250]; >> Y=[13.681,13.526,13.098,13.095]; >> [A,C,L,wcgs,Cw]=newpoly(X,Y) 输出3阶牛顿插值多项式L 及其系数向量C 差商的矩阵A ，插值余项wcgs 及其 ) ()()1(ξ+n n f x R 的系数向量Cw 。 A = 13.6810 0 0 0 13.5260 -0.1550 0 0 13.0980 -0.0713 0.0120 0 13.0950 -0.0030 0.0098 -0.0003 C = 1.0e+003 *
第06章_MATLAB数值计算_例题源程序汇总
第6章 MATLAB 数值计算 例6.1 求矩阵A 的每行及每列的最大和最小元素，并求整个矩阵的最大和最小元素。 1356 78256323578255631 01-???? -? ?=???? -??A A=[13,-56,78;25,63,-235;78,25,563;1,0,-1]; max(A,[],2) %求每行最大元素 min(A,[],2) %求每行最小元素 max(A) %求每列最大元素 min(A) %求每列最小元素 max(max(A)) %求整个矩阵的最大元素。也可使用命令：max(A(:)) min(min(A)) %求整个矩阵的最小元素。也可使用命令：min(A(:)) 例6.2 求矩阵A 的每行元素的乘积和全部元素的乘积。 A=[1,2,3,4;5,6,7,8;9,10,11,12]; S=prod(A,2) prod(S) %求A 的全部元素的乘积。也可以使用命令prod(A(:)) 例6.3 求向量X =(1！，2！，3！，…，10！)。 X=cumprod(1:10) 例6.4 对二维矩阵x ，从不同维方向求出其标准方差。 x=[4,5,6;1,4,8] %产生一个二维矩阵x y1=std(x,0,1) y2=std(x,1,1) y3=std(x,0,2) y4=std(x,1,2) 例6.5 生成满足正态分布的10000×5随机矩阵，然后求各列元素的均值和标准方差，再求这5列随机数据的相关系数矩阵。 X=randn(10000,5); M=mean(X) D=std(X) R=corrcoef(X)
例6.6 对下列矩阵做各种排序。 185412613713-?? ??=?? ??-?? A A=[1,-8,5;4,12,6;13,7,-13]; sort(A) %对A 的每列按升序排序 -sort(-A,2) %对A 的每行按降序排序 [X,I]=sort(A) %对A 按列排序，并将每个元素所在行号送矩阵I 例6.7 给出概率积分 2 (d x x f x x -? e 的数据表如表6.1所示，用不同的插值方法计算f (0.472)。 x=0.46:0.01:0.49; %给出x ，f(x) f=[0.4846555,0.4937542,0.5027498,0.5116683]; format long interp1(x,f,0.472) %用默认方法，即线性插值方法计算f(x) interp1(x,f,0.472,'nearest') %用最近点插值方法计算f(x) interp1(x,f,0.472,'spline') %用3次样条插值方法计算f(x) interp1(x,f,0.472,'cubic') %用3次多项式插值方法计算f(x) format short 例6.8 某检测参数f 随时间t 的采样结果如表6.2，用数据插值法计算t =2，7，12，17，22，17，32，37，42，47，52，57时的f 值。 T=0:5:65; X=2:5:57;
《有限元基础教程》_【MATLAB算例】3.3.7(2)__三梁平面框架结构的有限元分析(Beam2D2Node)
【MA TLAB 算例】3.3.7(2) 三梁平面框架结构的有限元分析 (Beam2D2Node) 如图3-19所示的框架结构，其顶端受均布力作用，结构中各个 截面的参数都为：113.010Pa E =?，746.510I m -=?，426.810A m -=?。试基 于MA TLAB 平台求解该结构的节点位移以及支反力。 图3-19 框架结构受一均布力作用 解答：对该问题进行有限元分析的过程如下。 （1） 结构的离散化与编号 将该结构离散为3个单元，节点位移及单元编号如图3-20所示， 有关节点和单元的信息见表3-5。 （a ） 节点位移及单元编号
（b）等效在节点上的外力 图3-20 单元划分、节点位移及节点上的外载 （2）各个单元的描述 首先在MA TLAB环境下，输入弹性模量E、横截面积A、惯性矩I和长度L，然后针对单元1，单元2和单元3，分别二次调用函数Beam2D2Node_ElementStiffness，就可以得到单元的刚度矩阵k1(6×6)和k2(6×6)，且单元2和单元3的刚度矩阵相同。 >> E=3E11; >> I=6.5E-7; >> A=6.8E-4; >> L1=1.44; >> L2=0.96; >> k1=Beam2D2Node_Stiffness(E,I,A,L1); >> k2=Beam2D2Node_Stiffness(E,I,A,L2); （3）建立整体刚度方程 将单元2和单元3的刚度矩阵转换成整体坐标下的形式。由于该结构共有4个节点，则总共的自由度数为12，因此，结构总的刚度矩阵为KK(12×12)，对KK清零，然后两次调用函数Beam2D2Node_Assemble进行刚度矩阵的组装。 >> T=[0,1,0,0,0,0;-1,0,0,0,0,0;0,0,1,0,0,0;0,0,0,0,1,0;0,0,0,-1,0,0;0,0,0,0,0,1] ; >> k3=T'*k2*T; >> KK=zeros(12,12); >> KK=Beam2D2Node_Assemble(KK,k1,1,2);
数值分析的MATLAB程序
列主元法 function lianzhuyuan(A,b) n=input('请输入n:') %选择阶数A=zeros(n,n); %系数矩阵A b=zeros(n,1); %矩阵b X=zeros(n,1); %解X for i=1:n for j=1:n A(i,j)=(1/(i+j-1)); %生成hilbert矩阵A end b(i,1)=sum(A(i,:)); %生成矩阵b end for i=1:n-1 j=i; top=max(abs(A(i:n,j))); %列主元 k=j; while abs(A(k,j))~=top %列主元所在行 k=k+1; end for z=1:n %交换主元所在行a1=A(i,z); A(i,z)=A(k,z); A(k,z)=a1; end a2=b(i,1); b(i,1)=b(k,1); b(k,1)=a2; for s=i+1:n %消去算法开始m=A(s,j)/A(i,j); %化简为上三角矩阵 A(s,j)=0; for p=i+1:n A(s,p)=A(s,p)-m*A(i,p); end b(s,1)=b(s,1)-m*b(i,1); end end X(n,1)=b(n,1)/A(n,n); %回代开始 for i=n-1:-1:1 s=0; %初始化s for j=i+1:n s=s+A(i,j)*X(j,1);
end X(i,1)=(b(i,1)-s)/A(i,i); end X 欧拉法 clc clear % 欧拉法 p=10; %贝塔的取值 T=10; %t取值的上限 y1=1; %y1的初值 r1=1; %y2的初值 %输入步长h的值 h=input('欧拉法please input number(h=1 0.5 0.25 0.125 0.0625):h=') ; if h>1 or h<0 break end S1=0:T/h; S2=0:T/h; S3=0:T/h; S4=0:T/h; i=1; % 迭代过程 for t=0:h:T Y=(exp(-t)); R=(1/(p-1))*exp(-t)+((p-2)/(p-1))*exp(-p*t); y=y1+h*(-y1); y1=y; r=r1+h*(y1-p*r1); r1=r; S1(i)=Y; S2(i)=R; S3(i)=y; S4(i)=r; i=i+1; end t=[0:h:T]; % 红线为解析解，'x'为数值解 plot(t,S1,'r',t,S3,'x')
matlab实现数值分析插值及积分
Matlab实现数值分析插值及积分 摘要： 数值分析(numerical analysis)是研究分析用计算机求解数学计算问题的数值计算方法及其理论的学科，是数学的一个分支，它以数字计算机求解数学问题的理论和方法为研究对象。在实际生产实践中，常常将实际问题转化为数学模型来解决，这个过程就是数学建模。学习数值分析这门课程可以让我们学到很多的数学建模方法。 分别运用matlab数学软件编程来解决插值问题和数值积分问题。题目中的要计算差值和积分，对于问题一，可以分别利用朗格朗日插值公式，牛顿插值公式，埃特金逐次线性插值公式来进行编程求解，具体matlab代码见正文。编程求解出来的结果为：=+。 其中Aitken插值计算的结果图如下： 对于问题二，可以分别利用复化梯形公式，复化的辛卜生公式，复化的柯特斯公式编写程序来进行求解，具体matlab代码见正文。编程求解出来的结果为： 0.6932 其中复化梯形公式计算的结果图如下：
问题重述 问题一：已知列表函数 表格 1 分别用拉格朗日，牛顿，埃特金插值方法计算。 问题二：用复化的梯形公式，复化的辛卜生公式，复化的柯特斯公式计算积分，使精度小于5。 问题解决 问题一：插值方法 对于问题一，用三种差值方法：拉格朗日，牛顿，埃特金差值方法来解决。 一、拉格朗日插值法： 拉格朗日插值多项式如下： 首先构造1+n 个插值节点n x x x ,,,10 上的n 插值基函数，对任一点i x 所对应的插值基函数 )(x l i ，由于在所有),,1,1,,1,0(n i i j x j +-=取零值，因此)(x l i 有因子 )())(()(110n i i x x x x x x x x ----+- 。又因)(x l i 是一个次数不超过n 的多项式，所以只 可能相差一个常数因子，固)(x l i 可表示成： )())(()()(110n i i i x x x x x x x x A x l ----=+- 利用1)(=i i x l 得：
同济大学数值分析matlab编程题汇编
MATLAB 编程题库 1.下面的数据表近似地满足函数2 1cx b ax y ++=，请适当变换成为线性最小二乘问题，编程求最好的系数c b a ,,，并在同一个图上画出所有数据和函数图像. 625 .0718.0801.0823.0802.0687.0606.0356.0995 .0628.0544.0008.0213.0362.0586.0931.0i i y x ---- 解： x=[-0.931 -0.586 -0.362 -0.213 0.008 0.544 0.628 0.995]'; y=[0.356 0.606 0.687 0.802 0.823 0.801 0.718 0.625]'; A=[x ones(8,1) -x.^2.*y]; z=A\y; a=z(1); b=z(2); c=z(3); xh=-1:0.1:1; yh=(a.*xh+b)./(1+c.*xh.^2); plot(x,y,'r+',xh,yh,'b*')
2.若在Matlab工作目录下已经有如下两个函数文件，写一个割线法程序，求出这两个函数 10 的近似根，并写出调用方式： 精度为10 解： >> edit gexianfa.m function [x iter]=gexianfa(f,x0,x1,tol) iter=0; while(norm(x1-x0)>tol) iter=iter+1; x=x1-feval(f,x1).*(x1-x0)./(feval(f,x1)-feval(f,x0)); x0=x1;x1=x; end >> edit f.m function v=f(x) v=x.*log(x)-1; >> edit g.m function z=g(y) z=y.^5+y-1; >> [x1 iter1]=gexianfa('f',1,3,1e-10) x1 = 1.7632 iter1 = 6 >> [x2 iter2]=gexianfa('g',0,1,1e-10) x2 = 0.7549 iter2 = 8
(完整版)有限元大作业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 程序的原理如框图：
matlab数值分析例题
1、 在MATLAB 中用Jacobi 迭代法讨论线性方程组， 1231231234748212515 x x x x x x x x x -+=?? -+=-??-++=? （1）给出Jacobi 迭代法的迭代方程，并判定Jacobi 迭代法求解此方程组是否收敛。 （2）若收敛，编程求解该线性方程组。 解（1）：A=[4 -1 1;4 -8 1;-2 1 5] %线性方程组系数矩阵 A = 4 -1 1 4 -8 1 -2 1 5 >> D=diag(diag(A)) D = 4 0 0 0 -8 0 0 0 5 >> L=-tril(A,-1) % A 的下三角矩阵 L = 0 0 0 -4 0 0 2 -1 0 >> U=-triu(A,1) % A 的上三角矩阵 U = 0 1 -1 0 0 -1 0 0 0 B=inv(D)*(L+U) % B 为雅可比迭代矩阵 B = 0 0.2500 -0.2500 0.5000 0 0.1250 0.4000 -0.2000 0 >> r=eigs(B,1) %B 的谱半径
r = 0.3347 < 1 Jacobi迭代法收敛。 (2)在matlab上编写程序如下： A=[4 -1 1;4 -8 1;-2 1 5]; >> b=[7 -21 15]'; >> x0=[0 0 0]'; >> [x,k]=jacobi(A,b,x0,1e-7) x = 2.0000 4.0000 3.0000 k = 17 附jacobi迭代法的matlab程序如下： function [x,k]=jacobi(A,b,x0,eps) % 采用Jacobi迭代法求Ax=b的解 % A为系数矩阵 % b为常数向量 % x0为迭代初始向量 % eps为解的精度控制 max1= 300; %默认最多迭代300，超过300次给出警告D=diag(diag(A)); %求A的对角矩阵 L=-tril(A,-1); %求A的下三角阵 U=-triu(A,1); %求A的上三角阵 B=D\(L+U); f=D\b; x=B*x0+f; k=1; %迭代次数 while norm(x-x0)>=eps x0=x; x=B*x0+f; k=k+1; if(k>=max1) disp('迭代超过300次，方程组可能不收敛'); return; end end
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 - ?? ?? - ????-- ??？
数值分析matlab代码
1、%用牛顿法求f(x)=x-sin x 的零点，e=10^(-6) disp('牛顿法'); i=1; n0=180; p0=pi/3; tol=10^(-6); for i=1:n0 p=p0-(p0-sin(p0))/(1-cos(p0)); if abs(p-p0)<=10^(-6) disp('用牛顿法求得方程的根为') disp(p); disp('迭代次数为：') disp(i) break; end p0=p; end if i==n0&&~(abs(p-p0)<=10^(-6)) disp(n0) disp('次牛顿迭代后无法求出方程的解') end 2、disp('Steffensen加速'); p0=pi/3; for i=1:n0 p1=0.5*p0+0.5*cos(p0); p2=0.5*p1+0.5*cos(p1); p=p0-((p1-p0).^2)./(p2-2.*p1+p0); if abs(p-p0)<=10^(-6) disp('用Steffensen加速求得方程的根为') disp(p); disp('迭代次数为：') disp(i) break; end p0=p; end if i==n0&&~(abs(p-p0)<=10^(-6)) disp(n0) disp('次Steffensen加速后无法求出方程的解') end 1、%使用二分法找到方程 600 x^4 -550 x^3 +200 x^2 -20 x -1 =0 在区间[0.1,1]上的根， %误差限为 e=10^-4 disp('二分法')
a=0.2;b=0.26; tol=0.0001; n0=10; fa=600*(a.^4)-550*(a.^3)+200*(a.^2)-20*a-1; for i=1:n0 p=(a+b)/2; fp=600*(p.^4)-550*(p.^3)+200*(p.^2)-20*p-1; if fp==0||(abs((b-a)/2) roots([1 -1 -1]) x=linspace(0,2*pi,10); y=sin(x); xi=linspace(0,2*pi,100); y1=interp1(x,y,xi); y2=interp1(x,y,xi,'spline'); y3=interp1(x,y,xi,'cublic'); plot(x,y,'o',xi,y1,xi,y2,xi,y3) x=[0 300 600 1000 1500 2000]; y=[0.9689 0.9322 0.8969 0.8519 0.7989 0.7491]; xi=linspace(0,2000,20); yi=1.0332*exp(-(xi+500)/7756); y1=interp1(x,y,xi,'spline'); subplot(2,1,1);plot(x,y,'o',xi,yi,xi,y1,'*') p=polyfit(x,y,2); y2=polyval(p,xi); subplot(2,1,2);plot(x,y,'o',xi,yi,xi,y2,'*') x=[0 300 600 1000 1500 2000]; y=[0.9689 0.9322 0.8969 0.8519 0.7989 0.7491]; xi=linspace(0,2000,20); y1=interp1(x,y,xi,'spline'); subplot(2,1,1);plot(x,y,'-o', xi,y1,'-*') p=polyfit(x,y,2); y2=polyval(p,xi); subplot(2,1,2);plot(x,y,'-o',xi,y2,'-*') 姓名：刘刚学号：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) 2015.1.9 上机作业题报告 JONMMX 2000 1.Chapter 1 1.1题目 设S N =∑1j 2?1 N j=2 ，其精确值为 )1 1 123(21+--N N 。 （1）编制按从大到小的顺序1 1 131121222-+ ??+-+-=N S N ，计算S N 的通用程序。 （2）编制按从小到大的顺序1 21 1)1(111222-+ ??+--+-= N N S N ，计算S N 的通用程序。 （3）按两种顺序分别计算64210,10,10S S S ,并指出有效位数。（编制程序时用单精度） （4）通过本次上机题，你明白了什么？ 1.2程序 1.3运行结果 1.4结果分析 按从大到小的顺序，有效位数分别为：6，4，3。 按从小到大的顺序，有效位数分别为：5，6，6。 可以看出，不同的算法造成的误差限是不同的，好的算法可以让结果更加精确。当采用从大到小的顺序累加的算法时，误差限随着N 的增大而增大，可见在累加的过程中，误差在放大，造成结果的误差较大。因此，采取从小到大的顺序累加得到的结果更加精确。 2.Chapter 2 2.1题目 （1）给定初值0x 及容许误差ε，编制牛顿法解方程f(x)=0的通用程序。 （2）给定方程03 )(3 =-=x x x f ,易知其有三个根3,0,3321= *=*-=*x x x ○1由牛顿方法的局部收敛性可知存在,0>δ当),(0δδ+-∈x 时，Newton 迭代序列收敛于根x2*。试确定尽可能大的δ。 ○2试取若干初始值，观察当),1(),1,(),,(),,1(),1,(0+∞+-----∞∈δδδδx 时Newton 序列的收敛性以及收敛于哪一个根。 （3）通过本上机题，你明白了什么？ 2.2程序 数值分析幂法与反幂法 matlab程序 随机产生一对称矩阵，对不同的原点位移和初值(至少取3个)分别使用幂法求计算矩阵的主特征值及主特征向量，用反幂法求计算矩阵的按模最小特征值及特征向量。 要求 1）比较不同的原点位移和初值说明收敛性 2）给出迭代结果，生成DOC文件。 3）程序清单，生成M文件。 解答： >> A=rand(5) %随机产生5*5矩阵求随机矩阵 A = 0.7094 0.1626 0.5853 0.6991 0.1493 0.7547 0.1190 0.2238 0.8909 0.2575 0.2760 0.4984 0.7513 0.9593 0.8407 0.6797 0.9597 0.2551 0.5472 0.2543 0.6551 0.3404 0.5060 0.1386 0.8143 >> B=A+A' %A矩阵和A的转置相加，得到随机对称矩阵B B = 1.4187 0.9173 0.8613 1.3788 0.8044 0.9173 0.2380 0.7222 1.8506 0.5979 0.8613 0.7222 1.5025 1.2144 1.3467 1.3788 1.8506 1.2144 1.0944 0.3929 0.8044 0.5979 1.3467 0.3929 1.6286 B=?? ????? ???? ?? ???6286.13929.03467.15979.08044 .03929.00944 .12144.18506 .13788.13467.12144.15025.17222.08613.05979.08506.17222.02380.09173.08044.03788.18613 .09173 .04187.1 编写幂法、反幂法程序： function [m,u,index,k]=pow(A,u,ep,it_max) % 求矩阵最大特征值的幂法，其中 % A 为矩阵； % ep 为精度要求，缺省为1e-5； % it_max 为最大迭代次数，缺省为100； % m 为绝对值最大的特征值； % u 为对应最大特征值的特征向量； % index ，当index=1时，迭代成功，当index=0时，迭代失败 if nargin<4 it_max=100; end if nargin<3 ep=1e-5; end n=length(A); index=0; k=0; m1=0; m0=0.01; % 修改移位参数，原点移位法加速收敛，为0时，即为幂法 I=eye(n) T=A-m0*I while k<=it_max v=T*u; [vmax,i]=max(abs(v)); m=v(i); u=v/m; if abs(m-m1) 数值分析matlab实现高斯消元法： function[RA,RB,n,X]=gaus(A,b) B=[A b];n=length(b);RA=rank(A); RB=rank(B);zhica=RB-RA; if zhica>0, disp('请注意：因为RA~=RB，所以此方程组无解.') return end if RA==RB if RA==n disp('请注意：因为RA=RB=n，所以此方程组有唯一解.') X=zeros(n,1);C=zeros(1,n+1); for p=1:n-1 for k=p+1:n m=B(k,p)/B(p,p);B(k,p:n+1)=B(k,p:n+1)-m*B(p,p:n+1); end end b=B(1:n,n+1);A=B(1:n,1:n);X(n)=b(n)/A(n,n); for q=n-1:-1:1 X(q)=(b(q)-sum(A(q,q+1:n)*X(q+1:n)))/A(q,q); end else disp('请注意：因为RA=RB Matlab大数值计算题目 1、统计附件1中的数据，对其中的数据划分区间，从0到50，每 10个单位一个区间，分为5个区间，统计每个区间的数量，画出柱状图。 Matlab程序： clear;clc;close all Data=xlsread('数据.xls'); Q=0:10:50; n=length(Data); m=length(Q); T=zeros(size(Q)); for s=1:n for t=1:m-1 if Data(s)>Q(t)&Data(s) 2、统计附件2中第二列数据中1至100每个数字出现的总次数， 附件2中第三列为每出现第二列数字所对应的次数，最后画出柱状图。 Matlab程序： clear;clc;close all Data=load('WEIBOIDWITHCOMMENTS.txt'); DATA=Data(:,2); t=Data(:,3); % m=max(DATA); m=100; T=zeros(m,1); for i=1:m data=DATA; data(data~=ones(size(data))*i)=0; data(data~=0)=1; n=data.*t; N=sum(n); T(i)=N; end bar(T) 3、 找到矩阵迷宫的通路,矩阵第1行第1列为迷宫的入口，第8行 第8列为迷宫的出口。（0表示路，1表示墙） 000 00000011 11010000 01010010 11010010 11010010 00011010 01000011 11110?????????????????????????? Matlab 程序： 主程序： clear all clc maze=[0,0,0,0,0,0,0,0; 0,1,1,1,1,0,1,0; 0,0,0,0,1,0,1,0; 0,1,0,0,0,0,1,0; 0,1,0,1,1,0,1,0; 0,1,0,0,0,0,1,1; 0,1,0,0,1,0,0,0; 0,1,1,1,1,1,1,0]; total=0; maze(1,1)=3;第3章 MATLAB数值计算-习题 答案
基于matlab的有限元法分析平面应力应变问题刘刚
东南大学-数值分析上机题作业-MATLAB版
数值分析幂法与反幂法-matlab程序
数值分析算法在matlab中的实现
Matlab大数值计算题目