线性代数实践(教师班第10讲).ppt
《线性代数实践(教师班第10讲).ppt》由会员分享,可在线阅读,更多相关《线性代数实践(教师班第10讲).ppt(63页珍藏版)》请在三一文库上搜索。
1、第十章 后续课矩阵建模举例,10.1 多项式插值问题,例:试求三次插值多项式 ,使曲线通过以下4个点: (0,3),(1,0),(2,-1),(3,6) 。 解:这4个点的坐标应满足三次多项式函数,代入后有: 程序: A=1,0,0,0;1,1,1,1;1,2,4,8;1,3,9,27; B=3;0;-1;6;a=AB 得到a=3,-2,-2,1T,即,例10.1的图形, 绘图程序 ezplot(3-2*t-2*t2+t3) hold on,grid on plot(0:3,3,0,-1,6,x) line(1.5,1.5,0,6 ) axis(-1,4,-2,8) 求t=1.5处的插值函数值
2、 t1=1.5; p1=3-2*t1-2*t12+t13 plot(t1,p1,o),图10.1 例10.1的插值曲线,高阶的多项式插值,在一般情况下,当给出函数在n+1个点上的值时,就可以用n次多项式 对它进行插值。如果给出的点数(即方程数)大于n+1,方程组成为超定的,因而没有一个能满足方程组的解,得出的曲线将是以最小二乘意义下的误差靠近各点,于是插值就变为拟合。 插值也不一定是自变量的多项式,比如圆锥曲线方程 虽然它有6个系数,若用a除以此方程两端,得到的将是有5个待定系数的方程。如果给出x-y平面上的5个点,就可以列出5个线性方程来确定这5个系数。,10.2 坐标测量仪测定的拟合,比如
3、为了测量一个圆锥形截面的半径,可在x-y平面内测量其圆周上n个点的坐标(xi,yi) (i=1,n),然后拟合出此截面的方程。 对于每一组数据(xi,yi),代入圆锥曲线方程,移项可得: n个点就有n个方程。其结构相同,只是数据不同。可以把数据写成列向量,然后用元素群运算一次列出所有的n个方程。,例10.2 圆锥截面方程的拟合,设测量了圆周上7个点,其x,y坐标如下: x = -3.000 -2.000 -1.000 0 1.000 2.000 3.000 y= 3.03 3.90 4.35 4.50 4.40 4.02 3.26 试求出此圆锥截面的方程,并求其最大最小直径。 列出程序如下:
4、x=-3:3; % 把x,y赋值为列向量 y=3.03,3.90,4.35,4.50,4.40,4.02,3.26; A=x.*y, y.*y, x, y, ones(size(x); B=-x.2; % 列出系数矩阵A和B c=inv(A*A)*A*B, % 解超定方程得出c,程序运行结果,将程序运行生成的参数写出,应为: 即截面的方程为:,例10.2 曲线,ezplot(x2+0.005*x*y+0.9214*y2-0.2228*x-0.4050*y-16.8537) 画圆锥曲线 hold on, plot(x,y,x) 画测量点 grid on, axis equal % x,y两方向取
5、同样比例尺 拟合的图形和给定的测量点, 如图所示。从工程角度看, 这样的测量点布局对拟合 精度不大有利。,估计圆直径的方法,设圆周方程为: c1,c2为圆心的坐标,r为半径。整理上述方程,得到 用n个测量点坐标(xi,yi)代入,得到 这是关于三个未知数的n个线性方程,所以是一个超定问题。解出就可得知这个最小二乘圆的圆心坐标和半径r的值:,例10.2a 曲线程序及运行结果,程序ag1002a x=-3:3; % 把x,y赋值为列向量 y=3.03,3.90,4.35,4.50,4.40,4.02,3.26; A=2*x, 2*y, ones(size(x) % 求出系数矩阵A,B B=x.2+
6、y.2 c=inv(A*A)*A*B, % 求超定方程的解,得出c r=sqrt(c(3)+c(1)2+c(2)2) % 由c求出r 程序运行的最后结果为: c = 0.1018 工件圆心的x坐标 0.4996 工件圆心的y坐标 15.7533 r = 4.0017 工件的半径r,按圆形估计的最小二乘图形,以上就是本题要解的超定线性方程组 快速绘制此圆形的语句为: ezplot(x-0.1018)2+(y-0.4996)2-4.00172=0),10.3 天体轨道测量拟合,在适当的极坐标中,天体的位置应满足下列方程: 其中为常数而是轨道的偏心率,对于椭圆,对于抛物线,而对于双曲线。对一个新发现
7、的天体的观测得到了表中的数据,试确定其轨道的性质,并预测此天体在弧度时的位置。 解:对任何一组测量数据i,ri,可列出以下对和两个变量的线性方程:,0,r,求和的程序,五个测量点共有五个方程,轨道只有两个未知数,这是一个超定问题。先用元素群运算一次列出5个方程的系数矩阵,再用超定方程解的公式。theta=0.88,1.1,1.42,1.77,2.14; 测定的theta 和r 数组r =3, 2.3, 1.65, 1.25, 1.01 A= ones(5,1), r.*cos(theta), B=r, 方程组系数矩阵 X=inv(A*A)*A*B % 解超定方程 rg=X(1)/(1-X(2)
8、*cos(0.46) % 求theta=0.46处的r 运行结果为: 即轨道方程为 ,rg=5.3098,程序ag1002a的数据和图形,% 画出极坐标中轨迹的语句, ezpolar(1.4509/(1-0.8111*cos(theta), hold on,polar(theta, r,x), axis(-4,8,-5,5) %画出测定点 得到图形如右。,10.4 静力学问题,静力学研究物体受力后的平衡方程。一个物体在平面上平衡,需要三个平衡方程。如果是N个物体相互作用下的平衡,那么方程的总数就是3N个。空间物体的平衡需要6个平衡方程。这个方程组通常都是线性的,所以求解就可以归结为矩阵方程。它
9、可以避免解单个方程和单个数,只要把系数矩阵输入程序,就轻而易举地同时得出所有的解。 例10.4 两杆系统受力图见图,设已知:G1=200; G2=100; L1= 2; L2=sqrt(2); theta1 =30*pi/180; theta2 = 45*pi/180; 求所示杆系的支撑反力Na,Nb,Nc。,例10.4的图和方程,对杆件1, X=0 Nax + Ncx = 0; Y=0 Nay + Ncy - G1 = 0; M=Ncy*L1*cos(theta1)-Ncx*L1*sin(theta1)-G1*L1/2*cos(theta1)=0; 对杆件2 X=0 Nbx - Ncx =
10、0; Y=0 Nby - Ncy - G2 = 0; M=Ncy*L2*cos(theta2)+Ncx*L2*sin(theta2)+G2*L2/2*cos(theta2)=0;,例10.4的矩阵模型,这是一组包含六个未知数Nax, Nay, Nbx, Nby, Ncx, Ncy的六个线性代数方程,写成矩阵方程AX+B,其中 A=1,0,0,0,1,0;0,1,0,0,0,1;0,0,0,0,-sin(theta1),. cos(theta1);0,0,1,0,-1,0;0,0,0,1,0,-1;0,0,0,0,. sin(theta2),cos(theta2); B=0; G1; G1/2*
11、cos(theta1); 0; G2; -G2/2*cos(theta2); X = AB % 求解线性方程组,程序ag1004的运行结果,这样求解的方法不仅适用于全部静力学题目,而且可用于材料力学和结构力学中的静力学和静不定问题。,10.5 材料力学的静不定问题,建模:A点因各杆变形而引起的x方向位移x,y方向位移y,由几何关系,得变形方程: 其中Ki为杆i的刚度系数。再加上两个力平衡方程: 共有n+2个方程,其中包含n个未知力和两个待求位移x和y,方程组适定可解.。,例10.5 三杆静不定桁架,图示三杆桁架,挂一重物P=3000,设L=2m,各杆的截面积分别为A1=200e-6, A2=3
12、00e-6, A3=400e-6, 材料的弹性模量E=200e9,求各杆受力的大小。 解:此时应有两个平衡方程和三个位移协调方程: -N1*cos(a1)-N2-N3*cos(a3)=0; N1*sin(a1)-N3*sin(a3)=P; N1/K1=dx*cos(a1)+dy*sin(a1) N2/K2 = dx N3/K3=dx*cos(a3)dy*sin(a3),程序ag1005编写,令X为列向量N1;N2;N3;dx;dy,把上述五个线性方程组列成D*X=B的矩阵形式,从而就可用MATLAB的左除语句XDB来求解。因此程序ag1005如下: % 设定参数,计算杆长,刚度系数等,分行给D
13、赋值 D(1,:) =cos(a1),cos(a2),cos(a3),0,0; D(2,:) =sin(a1),sin(a2),sin(a3),0,0; D(3,:) =1/K1,0,0,-cos(a1),-sin(a1); D(4,:) =0,1/K2,0,-cos(a2),-sin(a2); D(5,:) = 0,0,1/K3,-cos(a3),-sin(a3); B = P;0;0;0;0; % 给B赋值 format long, X = DB %解方程组,用长格式显示结果,程序ag1005运行结果,执行此程序,显示的结果为: X = 1763.40607065591(N1) 591.1
14、4251029634(N2) -2995.72429657297(N3) 0.00016949097(dx) 0.00001970475(dy) 若用普通格式显示,dy=0.0000,实际上它不是零,这可从N2不等于零推想。在MATLAB中一个矩阵中包含数值相差很大的元素时,就得采用format long来显示,或者单独显示各个元素X(1),X(5)。读者还可改变几根杆的刚度系数,看它们如何影响各杆力的分布。,10.6 二自由度机械振动,右图表示了由两个质量、弹簧及阻尼器构成的二自由度振动系统,今要在给定两个质量的初始位置和初始速度的情况下求系统的运动。,解:设x1和x2分别表示两个质量关于它
15、们平衡位置的偏差值,则此二自由度振动系统的一般方程为:,建模的矩阵形式,可写成矩阵形式: (2) 其中 (3) 这是一个四阶常微分方程组。给出它的初始条件(初始位置X0和初始速度Xd0): (4) 可以求出它的解,但用解析方法求解非常麻烦。如果C=0,即无阻尼情况,则系统可解耦为两种独立的振动模态,通常的书上只给出解耦简化后的解。,线性微分方程的矩阵解,基本思路是把原始方程化成四个一阶方程矩阵方程组: (5) 这个方程在初始条件Y=Y0下的解为 。用MATLAB表示为Y=Y0*expm(A*t)。其中expm表示把(A*t)看成矩阵来求其指数。只要把Y,Y0和A找到就行。 先把方程(2)写成两
16、个一阶矩阵方程: (6) 于是, (7) 因 故Y和Y0都是41单列变量;,例10.6 二维振动系统的数值例,设m1=1; m2=9; k1 = 4; k2=2; c1和c2可选。求在初始条件x0 = 1;0和 xd0 = 0;-1下,系统的输出x1,x2曲线。 根据上面的模型可以写出程序ag1006如下。 m1=1; m2=9; k1 = 4; k2=2; % 输入各原始参数 c1=input(c1=);c2=input(c2=); % 输入阻尼系数 x0 = 1;0; xd0 = 0;-1; tf= 50; dt=0.1; % 给出初始条件及时间向量 M = m1,0;0,m2; K =
17、k1+k2, -k2; -k2, k2; % 构成二阶参数矩阵 C = c1+c2, -c2; -c2, c2; A=zeros(2,2),eye(2);-MK,-MC; % 构成四阶参数矩阵 y0=x0;xd0; % 四元变量的初始条件 for i=1:round(tf/dt)+1 % 设定计算点,作循环计算 t(i)=dt*(i-1); y(:,i)=expm(A*t(i)*y0; % 循环计算矩阵指数 end,程序运行结果的曲线,运行此程序,输入c1=0.2,c2=0.5所得的结果见左图。从中可清楚地看到振动的两种模态。特别是x1的运动反映了两种模态的迭合。给出不同的初始条件,各模态的幅
18、度也会变化。输入c1=0,c2=0所得的结果见右图,这时两种振动模态分别出现在两个波形中,即已经简单地解耦了,振动模态分析,键入:p,lamda=eig(A)得到 特征根lamda由两组共轭复根组成,它的虚部是振动的角频率,它的实部是它的衰减系数。所以矩阵特征根反映了两种振动模态的特征,使我们更进一步理解了“特征”两字的物理意义。,10.6 交流稳态电路的复数矩阵,稳态交流电路方程是线性的,但其系数是复数,变量是平面上的向量(也是复数),所以得到的是复数线性方程组。虽然线性代数理论只证到实数,但MATLAB中有关线性代数的函数都对实数和复数同样有效,所以都可以用在复数条件下。,例10.7 如图
19、10-10所示电路,设Z1=-j250,Z2=250,Is=20A,负载ZL=500j 500,求负载电压。,例10.7的矩阵模型,MATLAB程序ag1007,Z1=-j*250;Z2=250;ki=0.5;Is=2+j*0; zL=500+j*500; % 设定元件参数 a11=1/Z1+1/Z2;a12=-1/Z2;a13=0; % 设定系数矩阵A a21=-1/Z2;a22=1/Z2-1/zL;a23=-ki; a31=1/Z1;a32=0;a33=-1; A=a11,a12,a13;a21,a22,a23;a31,a32,a33; B=1;0;0; % 设定系数矩阵B X=AB*Is
20、; % 求方程解X=Ua;Ub;I1 Ub=X(2) % 负载电压 absUb=abs(Ub) % 负载电压的幅度 angleUb=angle(Ub)*180/pi % 负载电压的相角(度) compass(Ub); % 画出Ub的向量图 % compass(Is;X) % 画出四个向量的图,程序运行结果为:,负载电压Ub=-2.5000e+002-7.5000e+002i 幅度absUb = 790.5694, 相角angleUb = -108.4349 Ub向量见图。 键入compass(Is;X),将 画出四个向量的图,因为Ua,Ub的数值比Is,I1大了上百倍, Is,I1在图中无法显
21、示,10.8 线性系统零输入响应的计算,线性时不变连续系统的特性可用常微分方程表示为: 求其零输入响应。 解:在零输入条件u=0时,等式右端为零。系统响应的通解为: 其中,p是特征方程的n个根组成的向量p1,p2,pn,其每个分量的系数Cn则由y及其各阶导数的初始条件y0, Dy0, , D(n-1)y0来确定。,代入初始条件得到的矩阵方程,初始条件数应该和常数数相等,由此构成一个确定C1,Cn的线性代数方程组,写成: 矩阵V由特征根向量p确定,这种矩阵称为范德蒙特矩阵。,求零输入响应程序ag1008,在MATLAB中,有生成它的函数vander(p)。它产生的矩阵与本题的矩阵排列转了90度,
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 线性代数 实践 教师 10
链接地址:https://www.31doc.com/p-2122658.html