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;

f7='(xc2-xd)^2+(yc2-yd)^2=(xc1-xd)^2+(yc1-yd)^2';

f8='(xc3-xd)^2+(yc3-yd)^2=(xc1-xd)^2+(yc1-yd)^2'; %前两个机构方程

f9='xc2=cos(sita12)*xc1-sin(sita12)*yc1+xp2-xp1*cos(sita12)+yp1*sin(sita12)';

f10='yc2=sin(sita12)*xc1+cos(sita12)*yc1+yp2-xp1*sin(sita12)-yp1*cos(sita12)'; %由第一个位移矩阵方程得出

f11='xc3=cos(sita13)*xc1-sin(sita13)*yc1+xp3-xp1*cos(sita13)+yp1*sin(sita13)';

f12='yc3=sin(sita13)*xc1+cos(sita13)*yc1+yp3-xp1*sin(sita13)-yp1*cos(sita13)'; %由第二个位移矩阵方程得出

f7=subs(f7,{xd,yd},{vxd,vyd});

f8=subs(f8,{xd,yd},{vxd,vyd});

f9=subs(f9,{xp1,xp2,yp1,sita12},{vxp1,vxp2,vyp1,vsita12});

f10=subs(f10,{xp1,yp1,yp2,sita12},{vxp1,vyp1,vyp2,vsita12});

f11=subs(f11,{xp1,xp3,yp1,sita13},{vxp1,vxp3,vyp1,vsita13});

f12=subs(f12,{xp1,yp1,yp3,sita13},{vxp1,vyp1,vyp3,vsita13}); %代入具体数值

[xc1,xc2,xc3,yc1,yc2,yc3]=solve(f7,f8,f9,f10,f11,f12); %解方程

vxc1=vpa(xc1);

vyc1=vpa(yc1);

vxc2=vpa(xc2);

vyc2=vpa(yc2);

vxc3=vpa(xc3);

vyc3=vpa(yc3);

(vxc1-vxd)^2+(vyc1-vyd)^2;

(vxc2-vxd)^2+(vyc2-vyd)^2;

(vxc3-vxd)^2+(vyc3-vyd)^2; %去掉这三行分号可验证C点三个位置是否距离D点相等

%最终答案

xb1,yb1,xc1,yc1

Lab=sqrt((vxb1-vxa)^2+(vyb1-vya)^2)

Lbc=sqrt((vxb1-vxc1)^2+(vyb1-vyc1)^2)

Lcd=sqrt((vxc1-vxd)^2+(vyc1-vyd)^2)

Lad=sqrt((vxa-vxd)^2+(vya-vyd)^2) %得到四杆长

'曲柄存在条件：'

%得出四杆长后计算得到

'可靠到位条件：'

[vxc1-vxb1,vyc1-vyb1]*[vxc1-vxd,vyc1-vyd]'

[vxc2-vxb2,vyc2-vyb2]*[vxc2-vxd,vyc2-vyd]'

[vxc3-vxb3,vyc3-vyb3]*[vxc3-vxd,vyc3-vyd]'

'顺序到位条件：'

%未完成

输出结果：

xb1 =

(-7-4*sin(2/9*pi)+4*cos(2/9*pi))/(4*cos(2/9*pi)+4*sin(2/9*pi)-5)

yb1 =

-1

xc1 =

-6*(27+24*sin(2/9*pi)-64*cos(2/9*pi))/(-31*cos(2/9*pi)-5+39*sin(2/9*pi))

yc1 =

-2*(72*cos(2/9*pi)-175+192*sin(2/9*pi))/(-31*cos(2/9*pi)-5+39*sin(2/9*pi))

Lab =

1.0288436025165976748172169832223

Lbc =

2.9872531417317691216303250912289

Lcd =

6.9831476545729886023199865357226

Lad =

5

ans =

曲柄存在条件：

ans =

可靠到位条件：

ans =

14.605219997928496422368168445525

ans =

19.799913716084881287517588922012

ans =

20.814756669957613005391246805307

ans =

顺序到位条件：

Matlab程序设计教程第二版刘卫国课后参考答案(供参考)

第二章 1 求下列表达式的值。 （1） w=sqrt(2)*(1+0.34245*10^(-6)) (2) a=3.5; b=5; c=-9.8; x=(2*pi*a+(b+c)/(pi+a*b*c)-exp(2))/tan(b+c)+a (3) a=3.32; b=-7.9; y=2*pi*a^(2)*[(1-pi/4)*b-(0.8333-pi/4)*a] (4) t=[2,1-3*i;5,-0.65]; z=1/2*exp(2*t)*log(t+sqrt(1+t^(2))) 2 求下列表达式 A=[-1,5,-4;0,7,8;3,61,7]; B=[8,3,-1;2,5,3;-3,2,0]; (1) A+6*B A^2-B+eye (2) A*B A.*B B.*A (3) A/B B\A (4) [A,B] [A([1,3],:);B^2] 3 根据已知，完成下列操作 （1） A=[23,10,-0.778,0;41,-45,65,5;32,5,0,32;6,-9.54,54,3.14]; K=find(A>10&A<25); A(K) （2） A=[23,10,-0.778,0;41,-45,65,5;32,5,0,32;6,-9.54,54,3.14]; B=A(1:3,:) C=A(:,1:2) D=A(2:4,3:4) E=B*C （3） E

2 用if语句 score=input('请输入成绩:'); if score>=90&&score<=100 disp('A'); elseif score>=80&&score<=89 disp('B'); elseif score>=70&&score<=79 disp('C'); elseif score>=60&&score<=69; disp('D'); elseif score<60&&score>=0; disp('E'); else disp('出错'); end 用switch语句 score=input('请输入成绩:'); switch fix(score/10) case {9,10} disp('A'); case {8} disp('B'); case {7} disp('C'); case {6} disp('D'); case {0,1,2,3,4,5} disp('E'); otherwise disp('出错'); end 第四章1题 1) X=0:10; Y=x-x.^3/6; P lot(x,y) 2) t=0:0.01:2*pi; x=8.*cos(t); y=4*sqrt(2).*sin(t); plot(x,y) 2题

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 解：平面桁架元是既有局部坐标又有总体坐标的二维有限元；对各结点

实验二 Matlab程序设计基本方法1

实验二Matlab程序设计基本方法 覃照乘自092 电气工程学院 一、实验目的： 1、熟悉MATLAB 程序编辑与设计环境 2、掌握各种编程语句语法规则及程序设计方法 3、函数文件的编写和设计 4、了解和熟悉跨空间变量传递和赋值 二、实验基本知识： ◆for循环结构 语法：for i＝初值：增量：终值 语句1 …… 语句n end 说明：1．i＝初值：终值，则增量为1。 2．初值、增量、终值可正可负，可以是整数，也可以是小数，只须符合数学逻辑。 ◆while 循环结构 语法：while 逻辑表达式 循环体语句 end 说明：1、whiIe结构依据逻辑表达式的值判断是否执行循环体语勾。若表达式的值为真，执行循环体语句一次、在反复执行时，每次都要进行判断。若表达 式的值为假，则程序执行end之后的语句。 2、为了避免因逻辑上的失误，而陷入死循环，建议在循环体语句的适当位 置加break语句、以便程序能正常执行。（执行循环体的次数不确定； 每一次执行循环体后，一定会改变while后面所跟关系式的值。） 3、while循环也可以嵌套、其结构如下：

while逻辑表达式1 循环体语句1 while逻辑表达式2 循环体语句2 end 循环体语句3 end ◆if-else-end分支结构 if 表达式1 语句1 else if 表达式2(可选) 语句2 else(可选) 语句3 end end 说明：1.if结构是一个条件分支语句，若满足表达式的条件，则往下执行；若不满足，则跳出if结构。 2．else if表达式2与else为可选项，这两条语句可依据具体情况取舍。 3．注意：每一个if都对应一个end，即有几个if，记就应有几个end。 ◆switch-case结构 语法：switch表达式 case常量表达式1 语句组1 case常量表达式2 语句组2 …… otherwise 语句组n end

有限元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求得。 总体刚度矩阵如下：

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软件丰富可靠的矩阵运算、数据处理、图形绘制等便利工具，可使得计算和图象一体化。对于结构

实验3MATLAB程序设计

1，编写M 函数实现求一个数是否为素数，再编写一主程序（脚本文件），要求通过键盘输入一个整数，然后调用判断素数函数，从而确定它是否素数。 x=input('请输入一个整数x:'); if myprime(x) disp('您输入的整数x是一个素数。') else disp('您输入的数x不是一个素数。') end function y=myprime(x) y=1; for i=2:fix(sqrt(x)) if mod(x,i)==0 y=0; end end 2，编写M 函数统计一数值中零的个数，然后编写脚本文件，实现统计从1—2007 中零的总个数。 function num=number0(a) %统计十进制数值中0的个数 sa=num2str(a);%将数值装化为字符串 num=length(find(sa=='0'));% ));%求取字符串中'0’的个数 y=0;

for a=1:2006 num=number0(a); y=num+y; end disp(y) 504 3，编写程序计算x∈[-3,3],字长0.01：并画出曲线x = -3:0.01:3; y=zeros(size(x)); for i = 1:length(x) if -3<= x(i)& x(i)<=-1 y(i)=(-x(i).^2-4*x(i)-3)/ 2; elseif -1<= x(i) & x(i)<=1 y(i)=-x(i).^2+1; elseif 1<=x(:,i)<=3 y(i)=(-x(i).^2+4*x(i)-3)/2; end end plot(x,y) -3-2-10123

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程序设计与应用课后实验六八九

实验六 高层绘图操作 %第一题： 程序代码如下： x=linspace(0,2*pi,101); y=(0.5+3*sin(x)./(1+x.^2)).*cos(x); plot(x,y) 01234567 -1 -0.5 0.5 1 1.5 %第二题： %（1） 程序代码如下： x=linspace(-2*pi,2*pi,100); y1=x.^2; y2=cos(2*x); y3=y1.*y2; plot(x,y1,'b-',x,y2,'r:',x,y3,'y--'); text(4,16,'\leftarrow y1=x^2'); text(6*pi/4,-1,'\downarrow y2=cos(2*x)'); text(-1.5*pi,-2.25*pi*pi,'\uparrow y3=y1*y2');

-8 -6 -4 -2 2 4 6 8 -30-20 -10 10 20 30 40 %（2） 程序代码如下： x=linspace(-2*pi,2*pi,100); y1=x.^2; y2=cos(2*x); y3=y1.*y2; subplot(1,3,1);%分区 plot(x,y1); title('y1=x^2');%设置标题 subplot(1,3,2); plot(x,y2); title('y2=cos(2*x)'); subplot(1,3,3); plot(x,y3); title('y3=x^2*cos(2*x)');

-10 10 0510 15202530 35 40y1=x 2 -10 10 -1-0.8 -0.6 -0.4-0.200.20.4 0.6 0.8 1y2=cos(2*x) -10 10 -30-20 -10 10 20 30 40 y3=x 2*cos(2*x) %（3） 程序代码如下： x=linspace(-2*pi,2*pi,20); y1=x.^2; subplot(2,2,1);%分区 bar(x,y1); title('y1=x^2的条形图');%设置标题 subplot(2,2,2); stairs(x,y1); title('y1=x^2的阶梯图'); subplot(2,2,3); stem(x,y1); title('y1=x^2的杆图'); subplot(2,2,4); fill(x,y1,'r');%如果少了'r'则会出错 title('y1=x^2的填充图'); %其他的函数照样做。

MATLAB程序设计实验报告

MATLAB 程序设计实验报告 一、实验目的 1. 通过实验熟悉MATLAB 仿真软件的使用方法； 2. 掌握用MATLAB 对连续信号时域分析、频域分析和s 域分析的方法，利用绘图命令绘制出典型信号的波形，了解这些信号的基本特征； 3. 掌握用MATLAB 对离散信号时域分析、频域分析和z 域分析的方法，利用绘图命令绘制出典型信号的波形，了解这些信号的基本特征； 4. 通过绘制信号运算结果的波形，了解这些信号运算对信号所起的作用。 二、实验设备 1. 计算机 ： 2. MATLAB R2007a 仿真软件 三、实验原理 对系统的时域分析 信号的时域运算包括信号的相加、相乘，信号的时域变换包括信号的平移、反折、倒相及信号的尺度变换。 （1）信号的相加和相乘：已知信号)(1t f 和)(2t f ，信号相加和相乘记为 )()(1t f t f =)(2t f +；)()(1 t f t f =)(2t f *。 （2）信号的微分和积分：对于连续时间信号，其微分运算是用diff 函数来完成的，其语句格式为：diff(function,’variable’,n)，其中function 表示需要进行求导运算的信号，或者被赋值的符号表达式；variable 为求导运算的独立变量；n 为求导的阶数，默认值为求一阶导数。连续信号的积分运算用int 函数来完成，语句格式为：diff(function,’variable’,a,b)，其中function 表示需要进行被积信号，或者被赋值的符号表达式；variable 为求导运算的独立变量；a,b 为积分上、下限，a 和b 省略时为求不定积分。 （3）信号的平移、翻转和尺度变换 信号的平移包含信号的左移与右移，信号的翻转包含信号的倒相与折叠，平移和翻转信号不会改变信号)(t f 的面积和能量。信号的尺度变换是对信号)(t f 在时间轴上的变化，可使信号压缩或扩展。)(at f 将原波形压缩a 倍，)/(a t f 将原波形扩大a 倍。 ￥ 对系统频率特性的分析

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 程序设计 实验目的： 1、掌握建立和执行M 文件的方法； 2、掌握实现选择结构的方法； 3、掌握实现循环结构的方法。 实验内容： 1. 编写用 5 次多项式拟合函数y=sin(x), x [0, 2 ]的脚本M 文件，要求绘图观察拟合的效果。 function shiyan1 x=0:0.5:2*pi y=sin(x) p=polyfit(x,y,5) x1=0:0.2:2*pi y1=polyval(p,x1) plot(x,y, 'b' ,x1,y1, '*r' x =

Columns 1 through 9 0 0.5000 1.0000 1.5000 2.0000 2.5000 3.0000 3.5000 4.0000 Columns 10 through 13 4.5000 5.0000 5.5000 6.0000 y = Columns 1 through 9 0 0.4794 0.8415 0.9975 0.9093 0.5985 0.1411 -0.3508 -0.7568 Columns 10 through 13 -0.9775 -0.9589 -0.7055 -0.2794 p = -0.0056 0.0881 -0.3967 0.2671 0.8902 0.0029 x1 = Columns 1 through 10 0 0.2000 0.4000 0.6000 0.8000 1.0000 1.2000 1.4000 1.6000 1.8000 Columns 11 through 20

2. 2.2000 2.4000 2.6000 2.8000 3.0000 3.2000 3.4000 3.6000 1.8001 Columns 21 through 30 4.0 4.2000 4.4000 4.6000 4.8000 5.0000 5.2000 5.4000 5.6000 5.8000 Columns 31 through 32 6.0 6.2000 y1 = Columns 1 through 10 0.29 0.1886 0.3786 0.5585 0.7172 0.8461 0.9391 0.9926 1.0048 0.9761 Columns 11 through 20 0.9083 0.8048 0.6701 0.5098 0.3301 0.1381 -0.0590 -0.2538 -0.4389 -0.6073 Columns 21 through 30 -0.7524 -0.8685 -0.9505 -0.9949 -0.9991 -0.9626 -0.8863 -0.7732 -0.6288 -0.4606 Columns 31 through 32

基于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 ??- ??? ???- ?? ? ???- ??? ?=??-?? ? ???---??? ???-??? ?

MATLAB程序设计及经典例题解析3

MATLAB程序设计 用MATLAB语言编写的程序，称为M文件。M文件可以根据调用方式的不同分为两类：命令文件(Script File)和函数文件(Function File)。 例3-1 分别建立命令文件和函数文件，将华氏温度f转换为摄氏温度c。 程序1：首先建立命令文件并以文件名f2c.m存盘。 clear; %清除工作空间中的变量 f=input('Input Fahrenheit temperature：'); c=5*(f-32)/9 然后在MATLAB的命令窗口中输入f2c，将会执行该命令文件，执行情况为： Input Fahrenheit temperature：73 c =22.7778 程序2：首先建立函数文件f2c.m。 function c=f2c(f) c=5*(f-32)/9 然后在MATLAB的命令窗口调用该函数文件。 clear; y=input('Input Fahrenheit temperature：'); x=f2c(y) 输出情况为： Input Fahrenheit temperature：70 c =21.1111 x =21.1111 3.1.2 M文件的建立与打开 M文件是一个文本文件，它可以用任何编辑程序来建立和编辑，而一般常用且最为方便的是使用MATLAB提供的文本编辑器。

1．建立新的M文件 为建立新的M文件，启动MATLAB文本编辑器有3种方法： (1) 菜单操作。从MATLAB主窗口的File菜单中选择New菜单项，再选择M-file命令，屏幕上将出现MATLAB 文本编辑器窗口。 (2) 命令操作。在MATLAB命令窗口输入命令edit，启动MATLAB文本编辑器后，输入M文件的内容并存盘。 (3) 命令按钮操作。单击MATLAB主窗口工具栏上的New M-File命令按钮，启动MATLAB文本编辑器后，输入M文件的内容并存盘。 2．打开已有的M文件 打开已有的M文件，也有3种方法： (1) 菜单操作。从MATLAB主窗口的File菜单中选择Open命令，则屏幕出现Open对话框，在Open对话框中选中所需打开的M文件。在文档窗口可以对打开的M文件进行编辑修改，编辑完成后，将M文件存盘。 (2) 命令操作。在MATLAB命令窗口输入命令：edit 文件名，则打开指定的M文件。 (3) 命令按钮操作。单击MATLAB主窗口工具栏上的Open File命令按钮，再从弹出的对话框中选择所需打开的M文件。 3.2 程序控制结构 3.2.1 顺序结构 1．数据的输入 从键盘输入数据，则可以使用input函数来进行，该函数的调用格式为： A=input(提示信息，选项)； 其中提示信息为一个字符串，用于提示用户输入什么样的数据。 如果在input函数调用时采用's'选项，则允许用户输入一个字符串。例如，想输入一个人的姓名，可采用命令： xm=input('What''s your name?','s'); 2．数据的输出 MATLAB提供的命令窗口输出函数主要有disp函数，其调用格式为

D实验五 M文件和MATLAB程序设计

实验五M文件和MATLAB程序设计 一、实验目的 matlab作为一种高级计算机语言，不仅可以命令行方式完成操作，也具有数据结构、控制流、输入输出等能力，本次实验通过熟悉和掌握m文件的建立与使用方法，以及函数与控制程序流程语句的使用，使学生具备一定的编程和程序调试能力。 1．掌握M文件的使用方法。 2．掌握if语句和switch语句的使用 3. 掌握循环语句的使用 4. 通过练习理解MATLAB编程方法。 二、实验原理 1．m文件 用matlab语言编写的程序，称为m文件。M文件根据调用方式的不同分为两类，命令文件（Script file）和函数文件（Function file）。区别? 2．程序控制结构 1）顺序结构 2）选择结构 （1）if语句a) 单分支if语句b) 双分支if语句c) 多分支if语句 （2）switch 语句 （3）try语句 3）循环结构 （1）for 语句 （2）while语句 （3）break语句、continue语句、return使用,区别? 3．函数文件 function 输出形参表＝函数名（输入形参表） 注释说明部分 函数体语句 注意事项？ 三、实验要求 1．首先上机练习PPT中各种流程控制语句的有关实例。 2．然后上机练习下面的实验习题。 四、实验习题

1．数论中一个有趣的题目：任意一个正整数，若为偶数，则用2除之，若为奇数，则与3相乘再加上1。重复此过程，最终得到的结果为1。如： 2→1 3→10→5→16→8→4→2→1 6→3→10→5→16→8→4→2→1 运行下面的程序，按程序提示输入n=1,2,3,5,7，8,9等数来验证这一结论。 %classic "3n+1" problem from number theory. while 1 n=input('Enter n,negative quits:'); if n<=0 break end a=n; while n>1 if rem(n,2)==0 n=n/2; else n=3*n+1; end a=[a,n]; end a end Enter n,negative quits: 2． 编程求满足∑=>m i i 11000020的最小m 值。 a=0; i=1; while (a<100000) a=a+pow2(i); i=i+1; end m=i-1 3． 编写一个函数，计算下面函数的值，给出x 的值，调用该函数后，返回y 的值。 function [y]=myfun1(x)

《结构力学习题集》-矩阵位移法习题及答案

第八章 矩阵位移法 – 老八校 一、判断题： 1、单元刚度矩阵反映了该单元杆端位移与杆端力之间的关系。 2、单元刚度矩阵均具有对称性和奇异性。 3、局部坐标系与整体坐标系之间的坐标变换矩阵T 是正交矩阵。 4、结构刚度矩阵反映了结构结点位移与荷载之间的关系。 5、结构刚度方程矩阵形式为：[]{}{}K P ?=，它是整个结构所应满足的变形条件。 6、图示结构用矩阵位移法计算时（计轴向变形）未知量数目为8个。 7、在直接刚度法的先处理法中，定位向量的物理意义是变形连续条件和位移边界条件。 8、等效结点荷载数值等于汇交于该结点所有固端力的代数和。 9、矩阵位移法中，等效结点荷载的“等效原则”是指与非结点荷载的结点位移相等。 10、矩阵位移法既能计算超静定结构，也能计算静定结构。 11、已知图示刚架各杆EI = 常数，当只考虑弯曲变形，且各杆单元类型相同时，采用先处理法进行结点位移编号，其正确编号是： (0,1,2) (0,0,0) (0,0,0) (0,1,3) (0,0,0)(1,2,0) (0,0,0)(0,0,3) (1,0,2) (0,0,0) (0,0,0)(1,0,3) (0,0,0) (0,1,2) (0,0,0)(0,3,4) A. B. C. D. 2134123412341234 ( )

二、计算题： 12、用先处理法计算图示结构刚度矩阵的元素133322,,K K K 。 12 3l l 4 l 5EI 2EI EA (0,0,0) (0,0,1) (0,2,3) (0,0,0) (0,2,4)(0,0,0) EI 13、用先处理法计算图示刚架结构刚度矩阵的元素153422,,K K K 。EI ，EA 均为常数。 l 14、计算图示结构整体刚度矩阵的元素665544,,K K K 。E 为常数。 l l 1 3 4 2 A , I A A /222A I , 2A 15、写出图示结构以子矩阵形式表达的结构原始刚度矩阵的子矩阵 [][]K K 22 24 ,。 [][]k k 1112 [][] k k 2122 [] k = i i i i i 单刚分块形式为 ：

矩阵位移法大作业

矩阵位移法大作业 学号：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程序设计实例

矩量法m ａtla ｂ程序设计实例： Ha ｌlen 方程求对称振子天线 一、条件与计算目标 已知: 对称振子天线长为Ｌ,半径为a ，且天线长度与波长得关系为,，设,半径a=0、０００0001，因此波数为。 目标： 用H ａll ｅn 方程算出半波振子、全波振子以及不同值得对应参数值。 求:（1）电流分布 （2）E 面方向图 (二维）,H 面方向图（二维）,半波振子空间方向性图(三维） 二、对称振子放置图 图1 半波振子得电流 分布 半波振子天线平行于z 轴放置，在ｘ轴与ｙ轴上得分量都为零，坐标选取方式有两种形式，一般选取图1得空间放置方 式。图１给出了天线得电流分布情况，由图可知，当天线很细时,电流分布近似正弦分布。 三、Ha ｌｌen 方程 得解题思路 ()()()()2 1 ' ' ' ' 12,cos sin sin 'z z i z z z z i z k z G z z dz c kz c kz E k z z dz j ωμ'++=-?? 对于中心馈电得偶极子，Hallen 方程为 ()22'1222 ('),'cos sin sin ,2L L i L L V i z G z z dz c kz c kz k z z j η + -- ++= <<+? 脉冲函数展开与点选配，得到 ()1121 ,''cos sin sin ,1,2,,2n n N z i n m m m m z n V I G z z dz c kz c kz k z m N j η +''=++= =???∑? 上式可以写成 矩阵形式为 四、结果与分析 （1）电流分布

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

结构静力计算与动力计算的对比分析 结构精力计算和结构动力计算是一个比较理论化和深度比较广的论述题目，在此，我仅凭本人有限的学识来展开对两者内容及关系的介绍和论述。也藉此契机，对结构力学上下册作一个比较系统的梳理和总结，为以后的学习以及工作打下坚实的基础。 首先，我想先介绍一下有关结构力学的基本概念，让读者可以带着一个整体、宏观的概念去深入理解具体的内部结构内容。 那么，我想从静力荷载和动力荷载的含义入题。静力荷载是指其大小、方向和位置不随时间变化或变化很缓慢的荷载，它不致使结构产生显著的加速度，因而可以略去惯性力的影响。结构的自重及其他恒荷载即属于静力荷载。动力荷载是指随时间迅速变化的荷载，它将引起结构振动，使结构产生不容忽视的及速度，因而必须考虑惯性力的影响。除荷载外，还有其他一些非荷载因素作用也可使结构产生内力和位移，例如温度变化、制造误差、材料收缩以及松弛、徐变等。 在结构静力计算中，最核心的内容就是计算结构的位移，而一切都要从虚功原理说起。虚功原理的两种表述：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 结 构