数值分析幂法与反幂法-matlab程序

数值分析幂法与反幂法-matlab程序
数值分析幂法与反幂法-matlab程序

数值分析幂法与反幂法 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)

m=m+m0; m1=m; k=k+1; end

function[m,u,index,k]=pow_inv(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;

% 修改移位参数,原点移位法加速收敛,为0时,即为反幂法I=eye(n);

T=A-m0*I;

invT=inv(T);

while k<=it_max

v=invT*u;

[vmax,i]=max(abs(v));

m=v(i);

u=v/m;

if abs(m-m1)

index=1;

break;

end

m1=m;

k=k+1;

end

m=1/m;

m=m+m0;

修改输入的m0的值,所得结果:

??5??

MATLAB在结构力学分析中的应用

MATLAB在结构力学分析中的应用 摘要:传统的手算方法解超静定结构工作量繁重,有时甚至是不可能,运用结构有限元编程的一般方法,通过两个实例的对照,展示MATLAB在结构力学分析中的应用,MATLAB具有高性能,方法具有普遍的适用性,实现弯矩图自动绘制。 关键词: MATLAB结构有限元弯矩图 Abstract:While using the traditional manual method to resolve complex statically indeterminate structures, it is heavy workloads, sometimes even impossible,using finite element programming of the general method, Based on two examples, This paper introduces a method of application of MATLAB in structure mechanics, MATLAB has the advantages of high performance, it can be applied to many kinds of structures, realization of automatic drawing bending moment diagram. Key words: MATLAB; Finite element; Bend moment diagram 引言 结构力学[3]中,常利用传统的力法与位移法求解超静定结构,力法是几何问题,位移法把复杂的几何图乘转化为代数运算,但它们基本未知量很多时,系数构成的矩阵计算巨大,两者都不能满足科研工作者的需要。应用MATLAB 软件丰富可靠的矩阵运算、数据处理、图形绘制等便利工具,可使得计算和图象一体化。对于结构力学计算是十分有利的工具。 1基本方法 MATLAB结构有限元编程的基本思路是先分后合,即将结构分成各个单元和节点,桁架与刚架已经离散化,对于连续系统这一步极其重要,然后进行单元分析,集成整体刚度矩阵,引入边界条件,最后解方程。在求解平面桁架结构,虽然结构简单,用手算可得各杆件的轴力,但重复的过程太多,现在使用MATLAB语言来编制有限元位移法的程序时,则编程的难度明显降低,对有限元位移法的概念的理解更加深入,编程所需时间也大大减少。 图1为一平面桁架,各杆E=70GPaA=0.004,试用矩阵位移法求解各杆轴力 图1 解:平面桁架元是既有局部坐标又有总体坐标的二维有限元;对各结点

幂法及反幂法

随机产生一对称矩阵,对不同的原点位移和初值(至少取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

北航数值分析大作业第一题幂法与反幂法

《数值分析》计算实习题目 第一题: 1. 算法设计方案 (1)1λ,501λ和s λ的值。 1)首先通过幂法求出按模最大的特征值λt1,然后根据λt1进行原点平移求出另一特征值λt2,比较两值大小,数值小的为所求最小特征值λ1,数值大的为是所求最大特征值λ501。 2)使用反幂法求λs ,其中需要解线性方程组。因为A 为带状线性方程组,此处采用LU 分解法解带状方程组。 (2)与140k λλμλ-5011=+k 最接近的特征值λik 。 通过带有原点平移的反幂法求出与数k μ最接近的特征值 λik 。 (3)2cond(A)和det A 。 1)1=n λλ2cond(A),其中1λ和n λ分别是按模最大和最小特征值。 2)利用步骤(1)中分解矩阵A 得出的LU 矩阵,L 为单位下三角阵,U 为上三角阵,其中U 矩阵的主对角线元素之积即为det A 。 由于A 的元素零元素较多,为节省储存量,将A 的元素存为6×501的数组中,程序中采用get_an_element()函数来从小数组中取出A 中的元素。 2.全部源程序 #include #include void init_a();//初始化A double get_an_element(int,int);//取A 中的元素函数 double powermethod(double);//原点平移的幂法 double inversepowermethod(double);//原点平移的反幂法 int presolve(double);//三角LU 分解 int solve(double [],double []);//解方程组 int max(int,int); int min(int,int); double (*u)[502]=new double[502][502];//上三角U 数组 double (*l)[502]=new double[502][502];//单位下三角L 数组 double a[6][502];//矩阵A int main() { int i,k; double lambdat1,lambdat2,lambda1,lambda501,lambdas,mu[40],det;

有限元MATLABWord版

MATLAB报告 Matlab程序求解简要过程如下: (1)求取单元节点位移提取矩阵T 单元节点位移提取矩阵T本质上是置换矩阵群中的一个,结果可将任意杂乱的节点顺序置换成统一的顺序。另一方面其作用是对单元刚度矩阵进行“升维操作”,将单元刚度矩阵统筹到整体刚度矩阵上来,便于对总体节点位移矩阵和支座反力进行求取。 本程序分析过程中对单元1的节点提取是按顺序编号1-2-3,对单元2的节点提取是按顺序编号2-3-4。单元1的节点位移提取矩阵如下:

单元2的节点位移提取矩阵如下: (2)求取单元几何矩阵B 单元1的节点按编号顺序1-2-3分别进行对几何函数矩阵或算子矩阵的bi逆时针操作,对ci顺时针操作;单元2的节点按编号顺序2-3-4分别进行对几何函数矩阵的bi顺时针操作,对ci逆时针操作.在MATLAB程序中通过mod()取模函数来达到对节点的顺时针或逆时针循环操作。 单元1的几何矩阵如下:

单元2的几何矩阵如下: (3)求取应力矩阵S 单元应力矩阵满足S=D*B,其中D为弹性矩阵,B为单元几何矩阵 各单元的弹性矩阵如下: 单元1的应力矩阵如下: 单元2的应力矩阵如下: (4)求取单元刚度矩阵K

单元刚度矩阵K满足公式K=B’*D*B*t*A,其中t为平面板的厚度,A为单元面积,且单元刚度矩阵为对称矩阵。 单元1的刚度矩阵如下: 单元2的刚度矩阵如下: (5)求取总体刚度矩阵sumKK 由上述步骤求得的单元刚度矩阵K利用单元虚功原理和刚度方程可导出K’*δ=f,其中δ为单元节点位移列阵,f为单元等效节点载荷列阵,为了能将各个单元刚度方程统一到一个整体,便需要步骤(1)的单元节点提取矩阵对单元刚度方程进行变换,将两个变换结果联立便得到总体刚度方程,其中也可得到总体刚度矩阵sumKK,且总体刚度矩阵可由sumKK=ΣT’*K*T求得。 总体刚度矩阵如下:

实验8 反幂法

《数值分析》实验8 一.实验名称:反幂法 二、实验目的: (1) 掌握求矩阵按模最小特征值及其对应的特征向量的规范化反幂法; (2) 掌握原点移位法。 三、实验要求 (1) 按照题目要求完成实验内容 (2) 写出相应的实验原理与C 语言程序 (3) 给出实验结果,结果分析 (4) 写出相应的实验报告 四、实验题目: 用反幂法求以下矩阵的指定特征值及其特征向量(迭代终止误差为1e-3): 41411014110?? ? ? ??? 接近2的特征值及其特征向量. 程序: #include #include int main(){ int n = 3, i, k, j, max_k = 1000; double y[max_k][n], x[max_k][n], max_x[max_k], z[n], u[n][n], l[n][n], t, err = 1e-5, norm_dx = 1,tt=12;//norm_dx:dx 的一范数,tt:为了求接近tt 的特征值, for (i = 0; i < n; i++) for (j = 0; j < n; j++){ u[i][j] = 0; (j == i) ? (l[i][j] = 1) : (l[i][j] = 0); } double a[][3] = { 4, 1, 4, 1, 10, 1, 4, 1, 10 }; for (i = 0; i < n; i++) { a[i][i]-=tt; x[0][i] = 1; } for (i = 0; i < n; i++){ for (k = i; k < n; k++){ t = 0; for (j = 0; j < i; j++) t += l[i][j] * u[j][k];

MATLAB在结构力学分析中应用

MATLAB在结构力学分析中的应用摘要:传统的手算方法解超静定结构工作量繁重,有时甚至是不可能,运用结构有限元编程的一般方法,通过两个实例的对照,展示matlab在结构力学分析中的应用,matlab具有高性能,方法具有普遍的适用性,实现弯矩图自动绘制。 关键词: matlab结构有限元弯矩图 abstract:while using the traditional manual method to resolve complex statically indeterminate structures, it is heavy workloads, sometimes even impossible,using finite element programming of the general method, based on two examples, this paper introduces a method of application of matlab in structure mechanics, matlab has the advantages of high performance, it can be applied to many kinds of structures, realization of automatic drawing bending moment diagram. key words: matlab; finite element; bend moment diagram 引言 结构力学[3]中,常利用传统的力法与位移法求解超静定结构,力法是几何问题,位移法把复杂的几何图乘转化为代数运算,但它们基本未知量很多时,系数构成的矩阵计算巨大,两者都不能满足科研工作者的需要。应用matlab软件丰富可靠的矩阵运算、数据处理、图形绘制等便利工具,可使得计算和图象一体化。对于结构

数值分析之幂法及反幂法C语言程序实例

数值分析之幂法及反幂法C 语言程序实例 1、算法设计方案: ①求1λ、501λ和s λ的值: s λ:s λ表示矩阵的按模最小特征值,为求得s λ直接对待求矩阵A 应用反幂法即可。 1λ、501λ:已知矩阵A 的特征值满足关系 1n λλ<< ,要求1λ、及501λ时,可 按如下方法求解: a . 对矩阵A 用幂法,求得按模最大的特征值1m λ。 b . 按平移量1m λ对矩阵A 进行原点平移得矩阵1m B A I λ=+,对矩阵B 用反幂法 求得B 的按模最小特征值2m λ。 c . 321m m m λλλ=- 则:113min(,)m m λλλ=,13max(,)n m m λλλ=即为所求。 ②求和A 的与数5011 140 k k λλμλ-=+最接近的特征值 ik λ(k=0,1,…39): 求矩阵A 的特征值中与k μ最接近的特征值的大小,采用原点平移的方法: 先求矩阵 B=A-k μI 对应的按模最小特征值k β,则k β+k μ即为矩阵A 与k μ最接近的特征值。 重复以上过程39次即可求得ik λ(k=0,1,…39)的值。 ③求A 的(谱范数)条件数2cond()A 和行列式det A : 在(1)中用反幂法求矩阵A 的按模最小特征值时,要用到Doolittle 分解方法,在Doolittle 分解完成后得到的两个矩阵分别为L 和U ,则A 的行列式可由U 阵求出,即:det(A)=det(U)。 求得det(A)不为0,因此A 为非奇异的实对称矩阵,则: max 2()s cond A λλ= ,max λ和s λ分别为模最大特征值与模最小特征值。

乘幂反幂法

当前位置:第7章>>第1节>>7.1.3 逆幂法 逆幂法是求实方阵按模最小特征值及相应的特征向量的一种反迭代方法。 1. 求A按模最小的特征值 设非奇异矩阵A的n 个特征值为,其相应的特征向量为e ,则的特征值为 其相应的特征向量仍为。 按模最大的特征值的倒数则为矩阵A按模最小的特征值。 利用乘幂法求按模最大的特征值。 任取初始非零初始向量,作迭代序列 它等价于(7.5) 我们可以通过反迭代过程,即解方程组. 求得. 当k 充分大时,则有 在实际计算中,为了减少运算量,先将矩阵A作三角分解A=LR 然后再求解方程组

2.求在附近的特征值 设与最接近的特征值为即有 作矩阵,它的特征值和相应的特征向量为 若用逆幂法于矩阵,则有 则可求出矩阵的按模最小的特征值和相应的特征向量为 于是得A在附近的特征值和相应的特征向量为 (7.6) 例3 用逆幂法求矩阵在3.4附近的特征值和相应的特征向量 解对进行三角分解得:

用半次迭代法,取,则 得 再解 得 再解 得 于是 练习7.1 1.用乘幂法求矩阵按模最大特征值与特征向量 . 乘幂法的计算公式

设矩阵A的n个特征值按模的大小排列为: 其相应的特征向量为且它们是线性无关的。 先任取非零初始向量,作迭代序列 首先将表示为 所以 为了得出计算和的公式,下面分三种情况讨论 1.为实根,且。 当不为0,k充分大时,则有 所以(7.1)2.为实根,且。 当不为0,k充分大时,则有

(7.2)于是得 从而有 (7.3)3.,且。当k充分大时,则有 在实际应用幂法时,可根据迭代向量个分量的变化情况判断属于那种情况。 若迭代向量各分量单调变化,且有关系式,则属于第1种情况; 若迭代向量各分量不是单调变化,但有关系式,则属于第2种情况;

MATLAB程序:已知三个位置设计平面四杆机构求解程序(位移矩阵法)

%MATLAB程序:已知三个位置设计平面四杆机构求解程序(位移矩阵法) clear;clc; %凡是变量名前带v的为数值变量,不带的是符号变量 vxp1=0; vyp1=0; vsita1=0*pi/180; vxp2=-2; vyp2=6; vsita2=40*pi/180; vxp3=-10; vyp3=8; vsita3=90*pi/180; %精确位置P1,P2,P3及各角度 vsita12=vsita2-vsita1; vsita13=vsita3-vsita1; vxa=-10; vya=-2; vxd=-5; vyd=-2; %选定A,D点 %所有数值均在此确定,更改此处即可解出不同数值的四杆机构位移矩阵方程 syms xp1 yp1 xp2 yp2 xp3 yp3 sita12 sita13; syms xa ya xb1 yb1 xb2 yb2 xb3 yb3; f1='(xb2-xa)^2+(yb2-ya)^2=(xb1-xa)^2+(yb1-ya)^2'; f2='(xb3-xa)^2+(yb3-ya)^2=(xb1-xa)^2+(yb1-ya)^2'; %前两个机构方程 f3='xb2=cos(sita12)*xb1-sin(sita12)*yb1+xp2-xp1*cos(sita12)+yp1*sin(sita12)'; f4='yb2=sin(sita12)*xb1+cos(sita12)*yb1+yp2-xp1*sin(sita12)-yp1*cos(sita12)'; %由第一个位移矩阵方程得出 f5='xb3=cos(sita13)*xb1-sin(sita13)*yb1+xp3-xp1*cos(sita13)+yp1*sin(sita13)'; f6='yb3=sin(sita13)*xb1+cos(sita13)*yb1+yp3-xp1*sin(sita13)-yp1*cos(sita13)'; %由第二个位移矩阵方程得出 f1=subs(f1,{xa,ya},{vxa,vya}); f2=subs(f2,{xa,ya},{vxa,vya}); f3=subs(f3,{xp1,xp2,yp1,sita12},{vxp1,vxp2,vyp1,vsita12}); f4=subs(f4,{xp1,yp1,yp2,sita12},{vxp1,vyp1,vyp2,vsita12}); f5=subs(f5,{xp1,xp3,yp1,sita13},{vxp1,vxp3,vyp1,vsita13}); f6=subs(f6,{xp1,yp1,yp3,sita13},{vxp1,vyp1,vyp3,vsita13}); %代入具体数值 [xb1,xb2,xb3,yb1,yb2,yb3]=solve(f1,f2,f3,f4,f5,f6); %解方程 vxb1=vpa(xb1); vyb1=vpa(yb1); vxb2=vpa(xb2); vyb2=vpa(yb2); vxb3=vpa(xb3); vyb3=vpa(yb3); (vxb1-vxa)^2+(vyb1-vya)^2; (vxb2-vxa)^2+(vyb2-vya)^2; (vxb3-vxa)^2+(vyb3-vya)^2; %去掉这三行分号可验证B点三个位置是否距离A点相等 syms xd yd xc1 yc1 xc2 yc2 xc3 yc3;

幂法和反幂法的matlab实现

幂法和反幂法的matlab实现

幂法求矩阵主特征值及对应特征向量 摘要 矩阵特征值的数值算法,在科学和工程技术中很多问题在数学上都归结为矩阵的特征值问题,所以说研究利用数学软件解决求特征值的问题是非常必要的。实际问题中,有时需要的并不是所有的特征根,而是最大最小的实特征根。称模最大的特征根为主特征值。 幂法是一种计算矩阵主特征值(矩阵按模最大的特征值)及对应特征向量的迭代方法,它最大的优点是方法简单,特别适用于大型稀疏矩阵,但有时收敛速度很慢。 用java来编写算法。这个程序主要分成了四个大部分:第一部分为将矩阵转化为线性方程组;第二部分为求特征向量的极大值;第三部分为求幂法函数块;第四部分为页面设计及事件处理。其基本流程为幂法函数块通过调用将矩阵转化为线性方程组的方法,再经过一系列的验证和迭代得到结果。

关键字:主特征值;特征向量;线性方程组;幂法函数块 POWER METHOD FOR FINDING THE EIGENVALUES AND CORRESPONDING EIGENVECTORS OF THE MATRIX ABSTRACT Numerical algorithm for the eigenvalue of matrix, in science and engineering technology, a

lot of problems in mathematics are attributed matrix characteristic value problem, so that studies using mathematical software to solve the eigenvalue problem is very necessary. In practical problems, sometimes need not all eigenvalues, but the maximum and minimum eigenvalue of real. The characteristic value of the largest eigenvalue of the modulus maximum. Power method is a calculation of main features of the matrix values (matrix according to the characteristics of the largest value) and the corresponding eigenvector of iterative method. It is the biggest advantage is simple method, especially for large sparse matrix, but sometimes the convergence speed is very slow. Using java to write algorithms. This program is divided into three parts: the first part is the matrix is transformed into linear equations; the second part for the sake of feature vector of the maximum; the third part is

数值分析幂法c语言实现

1.实验目的: 1熟练掌握C 语言程序设计,编程求解问题。 2.运用幂法求解住特征值和特征向量。 2.实验内容: 例题: 用幂法求 A= ??????????0.225.05.025.00.10.15.00.10.1 的特征值和特征向量。 完整代码以及截图如下: #include "stdio.h" #include "math.h" #define M 3 void main() { float fan(),max(),e1,e2,r1,r2; void au(),ex(),print_x(),std(); static float a[M][M]={{1.0,1.0,0.5},{1.0,1.0,0.25},{0.5,0.25,2.0}}; static float u0[M],u1[M],maxn0,maxn1; int i;

printf("*********************************\n"); printf("****** 幂法*********\n"); printf("******求特征值与特征向量*********\n"); printf("*********************************\n\n"); printf("input precision e1,e2:"); scanf("%f,%f",&e1,&e2); printf("\ninput u(%d):",M); for (i=0;ie1 || r2>e2) { printf("%4d",i++); print_x(u0); printf("\n"); ex(u0,u1); } else break; } while (1); } void au(a,u0,u1) float a[][M],u0[],u1[]; { int i,j; for (i=0;i

幂法反幂法求解矩阵大小特征值及其对应的特征向量

幂法反幂法求解矩阵大小特征值及其对应的特征向量

————————————————————————————————作者:————————————————————————————————日期:

数值计算解矩阵的按模最大最小特征值及对应的特征向量 一.幂法 1. 幂法简介: 当矩阵A 满足一定条件时,在工程中可用幂法计算其主特征值(按模最大)及其特征向量。矩阵A 需要满足的条件为: (1) 的特征值为A i n λλλλ,0||...||||21 ≥≥≥> (2) 存在n 个线性无关的特征向量,设为n x x x ,...,,21 1.1计算过程: i n i i i u x x αα,1 ) 0()0(∑==,有对任意向量不全为0,则有 1 11111221 12111 1 1 11 1 011)()(...u u a u a u λu λαu αA x A Ax x k n n k n k k n i i k i i n i i i k )(k (k))(k αλλλλλα++++=+=+++≈? ? ????+++======∑∑ 可见,当||1 2 λλ越小时,收敛越快;且当k 充分大时,有1)11 11)11111λαλαλ=??????==+++(k )(k k (k k )(k x x u x u x ,对应的特征向量即是)(k x 1+。 2 算法实现 . ,, 3,,1 , ).5() 5(,,,,||).4();max(,).3() (max(;0,1).2(,).1()() () (停机否则输出失败信息转置若转否则输出若计算最大迭代次数,误差限,初始向量输入矩阵βλβεβλβλε←+←<<-←←= ←←k k N k y x Ay x x abs x y k N x A k k k 3 matlab 程序代码

matlab连续梁程序的编制与使用

第三章连续梁程序的编制与使用 入结构力学领域中而产生的一种方法,而Matlab 语言正是进行矩阵运算的强大工具,因此,用 Matlab语言编写结构力学程序有更大的优越性。本 章将详细介绍如何利用Matlab语言编制连续梁结 构的计算程序。 矩阵位移法的解题思路是将结构离散为单元 (杆件),建立单元杆端力与杆端位移之间的关系 -单元刚度方程;再将各单元集成为原结构,在满 足变形连续条件和平衡条件时,建立整体刚度方 程;在边界条件处理完毕后,由整体刚度方程解出 节点位移,进而求出结构内力。 用矩阵位移法计算连续梁的步骤如下: 1)整理原始数据,如材料性质、荷载条件、约 束条件等,并进行编码:单元编码、结点编 码、结点位移编码、选取坐标系。 2)建立局部坐标系下的单元刚度矩阵。 3)建立整体坐标系下的单元刚度矩阵。 4)集成总刚。 5)建立整体结构的等效节点荷载和总荷载矩 阵 6)边界条件处理。 7)解方程,求出节点位移。 8)求出各单元的杆端内力。 实际上,上述步骤也是编制Matlab程序的基本 步骤,在求出计算结果后,还可以利用Matlab的绘图功能绘制结构图、内力图、 变形图等等。 图3-1程序流程图

3.1 程序说明 %******************************************************************* % 矩阵位移法解连续梁主程序 %******************************************************************* ●功能:运用矩阵位移法解连续梁的基本原理编制的计算主程序。 ●基本思想:结点(结点位移)编码默认为从左至右,从1开始顺序进行;杆 端弯矩的方向默认为逆时针。 ●荷载类型:可计算结点荷载,每单元作用的跨中集中力和均布荷载。 ●说明:主程序的作用是通过赋值语句、读取和写入文件、函数调用等完成算 法的全过程,即实现程序流程图的程序表达。 %----------------------------------------------------------------------------------------------------- 1 程序准备 format short e %设定输出类型 clear all %清除所有已定义变量 clc %清屏 ●说明: format short e -设定计算过程中显示在屏幕上的数字类型为短格式、科学计数法; clear all -清除所有已定义变量,目的是在本程序的运行过程中,不会发生变量名相同等可能使计算出错的情况; clc -清屏,使屏幕在本程序运行开始时 %----------------------------------------------------------------------------------------------------- 2 打开文件 FP1=fopen('input.txt','rt'); %打开输入数据文件存放初始数据FP2=fopen('output.txt','wt'); %打开输出数据文件存放计算结果 ●说明: FP1=fopen('input.txt','rt'); -打开已存在的输入数据文件input.txt,且设置其为只读格式,使程序在执行过程中不能改变输入文件中的数值,并用文件句柄FP1来 FP2=fopen('output.txt','wt'); -打开输出数据文件,该文件不存在时,通过此

数值分析幂法与反幂法-matlab程序

数值分析幂法与反幂法 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的平面刚架静力分析

基于MATLAB 的平面刚架静力分析 为了进一步理解有限元方法计算的过程,本文根据矩阵位移法的基本原理应用MATLAB 编制计算程序对以平面刚架结构进行了静力分析。本文还利用ANSYS 大型商用有限元分析软件对矩阵位移法的计算结果进行校核,发现两者计算结果相当吻合,验证了计算结果的可靠性。 一、 问题描述 如图1所示的平面刚架,各杆件的材料及截面均相同,E=210GPa ,截面为0.12×0.2m 的实心矩形,现要求解荷载作用下刚架的位移和内力。 5m 4m 3m 图1 二、矩阵位移法计算程序编制 为编制程序方便考虑,本文计算中采用“先处理法”。具体的计算步骤如下。

(1) 对结构进行离散化,对结点和单元进行编号,建立结构(整体)坐标系 和单元(局部)坐标系,并对结点位移进行编号; (2) 对结点位移分量进行编码,形成单元定位向量e λ; (3) 建立按结构整体编码顺序排列的结点位移列向量δ,计算固端力e F P 、等 效结点荷载E P 及综合结点荷载列向量P ; (4) 计算个单元局部坐标系的刚度矩阵,通过坐标变换矩阵T 形成整体坐标 系下的单元刚度矩阵e T e K T K T = ; (5) 利用单元定位向量形成结构刚度矩阵K ; (6) 按式1=K P δ- 求解未知结点位移; (7) 计算各单元的杆端力e F 。 根据上述步骤编制了平面刚架的分析程序。程序中单元刚度矩阵按下式计算。 32322 23 2 32 22 0000 1261260 064620 00001261260062640 EA EA l l EI EI EI EI l l l l EI EI EI EI l l l l K EA EA l l EI EI EI EI l l l l EI EI EI EI l l l l ??- ??? ???- ?? ? ???- ??? ?=??-?? ? ???---??? ???-??? ?

矩阵位移法大作业

矩阵位移法大作业 学号:151210122 姓名:谭逸天 班级:土木一班

编制原理: 使用Math Work公司开发的科学与工程计算机软件——MATLAB, 利用其矩阵运算的便利性,将题目要求结构的基本信息编入脚本命令文件中,并编入求解步骤。加上刚度信息的输入指令,以及提取解答要求信息并输出的指令。令使用者只需输入结构材料相关信息便可计算题目对应悬索—拱组合体系的信息,并直接在命令窗口输出。 利用计算套路的重复性,程序开发时进行模块化设计。再由重复单元完成多次、重复的运算。 从整体性考虑,数据储存采用“算后集装,装后回收”对变量及数组重复使用,由配音进行简单命名,提高可辨识度。由于计算套路及程序本身高度模块化,并且题目所需个体信息相对于整体极少,提取个体化的信息只需简单改造命令模块,从整体信息中提取处理得出。编程所需的“数据化”“编码”等预处理由人工在编程开始前完成,由左下斜索基座作原点,正右向为X轴正向,正上为Y轴正向,建立右手系。编码顺序从左倒右由上及下,并用先处理法处理基座。(如下图所示)

6 7 共45个单元,32个结点编号,71个位移编号。 本人学号对应节间数m=14;f1=7L/4;f2=7L/10;h=7L/2;以上数据 为编程中人工设定值,结构的其余信息根据用户的输入进行计算得出。

程序说明: 初始计算结构在坐标系中的坐标信息,手动编入悬索与拱的曲线关键点信息,代入方程求解。随后由循环语句模块计算并存储结构中各类杆件的角度、长度信息,采用以直代曲的方法处理曲线。 由于先处理法,两端各四个单元不与其余单元通用编码递进规律,采用单独的语句进行计算并集装入总体信息储存矩阵中,其余规律性单元信息由循环的语句模块进行集装,便于之后的计算。定位向量统一装至71行6列的矩阵“dingwei”中,单元的长度与夹角信息统一装至71行2列的矩阵“danyuan”中,第一列为长度,第二列为角度。使两个信息矩阵的行序号对应单元序号,便于之后使用。 之后进入单元分析部分。先是对上部悬索进行单元分析,此部分为桁架单元,从“danyuan”矩阵中提取长度信息与角度信息,结合 开始时输入的刚度信息组装单刚矩阵与坐标变换矩阵,进行坐标变换后直接提取定位向量进行集装部分总刚矩阵的步骤。集装命令通过循环嵌套配合判断语句,对单刚矩阵进行二维遍历,并提取合格的元素填充至对应位置。随后,通过少量改动实现对斜索、吊杆、拱、主塔的处理。 之后保留基本结构,进行单元结点荷载的分析,并集装出结构结点荷载矩阵。 之后通过简单矩阵运算即得结构结点位移列阵。 进入单元后处理。将集装循环语句进行改造,达成逆向提取单元结点位移的功能。提取之前存储的单元信息进行坐标变换。最后算出

数值方法课程设计幂法反幂法计算矩阵特征值和特征向量附Matlab程序

数值方法课程设计幂法反幂法计算矩阵特征值和特征向量附Matlab程序

矩阵的特征值与特征向量的计算 摘要 物理,力学,工程技术中的很多问题在数学上都归结于求矩阵特征值的问题,例如振动问题(桥梁的振动,机械的振动,电磁振动等)、物理学中某些临界值的确定问题以及理论物理中的一些问题。矩阵特征值的计算在矩阵计算中是一个很重要的部分,本文使用幂法和反幂法分别求矩阵的按模最大,按模最小特征向量及对应的特征值。 幂法是一种计算矩阵主特征值的一种迭代法,它最大的优点是方法简单,对于稀疏矩阵比较合适,但有时收敛速度很慢。其基本思想是任取一个非零的初始向量。由所求矩阵构造一向量序列。再经过所构造的向量序列求出特征值和特征向量。 反幂法用来计算矩阵按模最小特征向量及其特征值,及计算对应于一个给定近似特征值的特征向量。本文中主要使用反幂法计算一个矩阵的按模最小特征向量及其对应的特征值。计算矩阵按模最小特征向量的基本思想是将其转化为求逆矩阵的按模最大特征向量。然后经过这个按模最大的特征向量反推出原矩阵的按模最小特征向量。

关键词:矩阵;特征值;特征向量;冥法;反冥法 THE CALCULATIONS OF EIGENVALUE AND EIGENVECTOR OF MATRIX ABSTRACT Physics, mechanics, engineering technology in a lot of problems in mathematics are attributed to matrix eigenvalue problem, such as vibration (vibration of the bridge, mechanical vibration, electromagnetic vibration, etc.) in physics, some critical values determine problems and

结构静力计算和动力计算的对比分析

结构静力计算与动力计算的对比分析 结构精力计算和结构动力计算是一个比较理论化和深度比较广的论述题目,在此,我仅凭本人有限的学识来展开对两者内容及关系的介绍和论述。也藉此契机,对结构力学上下册作一个比较系统的梳理和总结,为以后的学习以及工作打下坚实的基础。 首先,我想先介绍一下有关结构力学的基本概念,让读者可以带着一个整体、宏观的概念去深入理解具体的内部结构内容。 那么,我想从静力荷载和动力荷载的含义入题。静力荷载是指其大小、方向和位置不随时间变化或变化很缓慢的荷载,它不致使结构产生显著的加速度,因而可以略去惯性力的影响。结构的自重及其他恒荷载即属于静力荷载。动力荷载是指随时间迅速变化的荷载,它将引起结构振动,使结构产生不容忽视的及速度,因而必须考虑惯性力的影响。除荷载外,还有其他一些非荷载因素作用也可使结构产生内力和位移,例如温度变化、制造误差、材料收缩以及松弛、徐变等。 在结构静力计算中,最核心的内容就是计算结构的位移,而一切都要从虚功原理说起。虚功原理的两种表述:1、对于刚体体系,刚体体系处于平衡的必要和充分条件是,对于任何虚位移,所有外力所作虚功总和为零;2、对于变形体系,其处于平衡的必要和充分条件是,对于任何虚位移,外力所作虚功总和等于各微段上的内力在其变形上所作的虚功总和,简单说,外力虚功等于内力虚功。 虚功方程: 由于力状态与位移状态是彼此独立无关的,因此运用单位荷载法: 由: 得位移计算一般公式: 同过几何关系可得弯矩图乘法便捷计算公式(为计算带来极大的方便): 力法: 力法典型方程: (系数δ?、的求解方法如同上述虚功原理的原理。) 该方程的物理意义为:基本结构在全部多余未知力和荷载共同作用下,在去掉各多余联系处沿各多余未知力方向的位移,应与原结构相应的位移相等。可见,力法可以求解出超静N u s s W F d Md F d ?γ=++∑∑∑???1k R N u s s F c F d Md F d ?γ?+=++∑∑∑∑??? N S S s k N s R F ds Md F d F M F F c EA EI GA γ?=++-∑∑∑∑???S w c Md A y M EI EI =∑?1111221211222200P P X X X X δδδδ++?=??++?=?基本体系 1 X 结 构

相关文档
最新文档