第三章 基本数值计算方法(一)

《第三章 基本数值计算方法(一)》由会员分享，可在线阅读，更多相关《第三章 基本数值计算方法(一)（32页珍藏版）》请在金锄头文库上搜索。

1、第三章 基本数值计算方法（一 ）1.线性代数的数值解法 2.矩阵特征值问题一、线性代数的数值解法线性代数解决的实际问题大体上就是三类:1）求线性代数方程组（包括欠定、适定 和超定）的解；2)分析向量的线性相关性；3)矩阵的特征值与对角化。 下面的将以求线性代数方程组为重点，围绕 这几个方面展开。一、线性代数的数值解法1. 线性代数方程组求解1.1 利用左除运算符的直接解法 对于线性方程组Ax=b，可以利用左除运算 符“”求解：x=Ab例3-1 用直接解法求解下列线性方程组。命令如下：A=2,1,-5,1;1,-5,0,7;0,2,1,-1;1,6,-1,-4; b=13,-9,6,0; x=A

2、b1.2 利用X=inv(A)*b语句的直接解法 方程AXb可改写为XA1 b 因此可用 X=inv(A)*b 解出。 例3-2A=2,1,-5,1;1,-5,0,7;0,2,1,-1;1,6,-1,-4;b=13,-9,6,0;X= X=inv(A)*b1.3最简行阶梯法解线性代数方程组问题的提出解：可以把线性方程组写成矩阵方程A*x=b，其中:A=6,1,6,-6;1,-1,9,9;-2,4,0,-4;4,2,7,-5; b=7;5;-7;-9 U=rref(A,b) 程序运行的结果为：这个矩阵代表了以下的等价方程组， 所以等式右端的四个数也就是方程的解。X1 = 15.58X2 = 14

3、.70X3 = -8.20X4 = 8.66 为什么要提出这种新的计算方法？ 把上例中第四个方程改为：，求其解。解：输入新参数 A=6,1,6,-6;1,-1,9,9;-2,4,0, -4;4,2,7,-778/222; b=7;5;-7;877/222; 键入U=rref(A,b)，得到这个最简行阶梯形式说明 原来的方程组是欠定的 。欠定方程组解的特点它等价于下列方程组： 这是一组包括四个变量的三个有效方程。因此没有唯一 的解。其中x4可以任意设定，即可以看作任意常数c， ：代入不同的c可以得到不同的解，因此欠定方程组有无数个解。这些解组成一根在空间中的直线。 用前两种方法试试？用方法一、二

4、：键入x=inv(A)*b或x=Ab，得到： Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 3.590822e -018.计算机告诉我们，这个结果是不准确的。其原因 在于矩阵A的行列式接近于零，不难检验， det(A)=0，rank(A)=3。就是说，A的逆极小。从 逆条件数RCOND = 3.590822e-018，可以得知计 算结果有效数位减小的程度。MATLAB的是有效 数位是16位，现在算得的结果有效数位要减去18 位，所以得出的结果是根本不能相信的。不相

5、容方程组的示例 如果把第四个方程的右端常数仍取为 -9，则其行阶梯变 换的结果为： 最后一个方程成了一个矛盾方程0=1。这说明方程组不 相容，无解。 由此也可以看出，线性方程组求解最好还是用行阶梯简 化的方法。因为它可以给出线性方程组的特征，避免计 算的盲目性。2.线性方程组应用举例2.1多项式插值问题 给出平面上4个点的坐标值如右表， (1)。求对它进行插值的三次多项式， (2)。求t=1.5处f的近似值。 (3)。如果要求此多项式多通过一点(-1,5)，求其系数。 解：用多项式 来插值，令它在四点上的值与表中相同，得到ti0123 f(ti )30-16多项式插值问题的解(1)这个矩阵方程

6、达到四阶，应该用计算机辅助求解了，编出 MATLAB程序exn554 C=1,0,0,0;1,1,1,1;1,2,4,8;1,3,9,27, b=3;0;-1;6,U=rref(C,b) 得到， 多项式为 多项式插值问题的解(2)（2）把t1=1. 5代入此多项式， 键入： f1=3-2*1.5-2*1.52+1.53 得到：f1=-1.125。 可以用以下语句画出插值图形 ， ezplot(t3-2*t2-2*t+3,- 1,4), grid on,hold on plot(0:3,3,0,-1,6,*) plot(1.5,-1.125,or) 可以得到图5-39，曲线通过 了图中四个给定的

7、插值点（ 用*号表示），圆圈为f(1.5) 的位置。多项式插值问题的解(3)(3)。若要此三次多项式多通过一点(-1,5)，将此点坐标代 入后得，方程组就成为五个方程四个未知数，很可能是 矛盾的，要靠计算来判断。用以下程序： A1=1,0,0,0;1,1,1,1;1,2,4,8;1,3,9,27;1,-1,1,-1, b1=3;0;-1;6;5,U01=rref(A1,b1) 得到 多项式插值问题的解(4)注意到系数增广矩阵最后一列是常数项，得出右方的同 解方程组。其最后一个方程成为0=1,说明方程组是不相 容的，属于超定方程，无通常意义下的解。 用MATLAB的pinv函数可以求出它的最小二

8、乘解，键入 ： a1=pinv(A1)*b1，得： 多项式成为多项式插值问题的解(5) 画出它的图形，可以继续键入 t=-1:0.1:4;plot(t,a1(1)+a1(2)*t+a1(3)*t.2+a1(4)*t.3,:r) plot(-1,5,x) 此曲线也用虚线画在图上。它不通过给定的5个点中的任 何一个，说明都有误差。误差向量E及五项误差的平方 和EE可以用矩阵方程计算如下 E=A1*a1-b1，EE=norm(E) 求得，结果为2.1化学反应方程的配平 配平下列小苏打与柠檬酸的反应。 解：建模，按每种物质的元素Na,H,C,O写成列向量， 按四种元素左右平衡列出四个方程，得：化学方程

9、配平程序 A=1,0,-3,0,0;1,8,0,-2,0;1,6,-6,0,-1;3,8,-7,-1,- 2,b=0;0;0;0 U=rref(A,b) 结果为 整数格式配平程序exn555运行结果 由此可见，x5必须最小取30，才能保证各个分 量xk为整数，于是得到 最后可以得出配平后的化学反应方程为 这个待配平的化学方程有四个方程和五个未知 数，所以它是一个欠定的线性方程组，其解乘 以任何常数仍然为方程的解。不过当我们另加 了一个条件系数取最小的正整数，结果就 是唯一的了。 二、求方阵的特征值和特征向量设有对称实矩阵，试求其特征方程、特征根和特征向量，并讨论矩阵的行 列式和迹与特征值的关系

10、。分三个步骤： 不用eig函数，按照原理，利用MATLAB函数分步求解； 用eig函数求解，进行对照校核； 求特征值矩阵的行列式与迹。例5-5-8的MATLAB程序exn558%（1）分步计算特征值和特征向量 A=2,4,9;4,2,4;9,4,18,% 输入矩阵参数 f=poly(A)% 求其特征多项式的系数向量f r=roots(f)% 求其特征根 for i=1:3 % 将三个特征根分别代入特征方程p(:,i)=null(A-r(i)*eye(3); % 求齐次方程的基础解 end, p,d=diag(r) % 列出特征向量矩阵和特征根矩阵%（2）用eig函数求特征值和特征向量： p1,

11、d1 = eig(A)% 一步求出特征值和特征向量 % 求原矩阵及特征根矩阵的行列式和迹 det(A), det(d), trace(A), trace(d), (1) 分步计算的结果(1) 分步计算的结果为： f = 1.0000 -22.0000 -37.0000 122.0000 r = 23.3603-3.06451.7042 p = 0.4149 -0.8483 0.32900.2419 0.4514 0.85890.8771 0.2767 -0.3925 说明此矩阵A的特征方程为，特征根为23.3603，- 3.0645，1.7042。将特征根在主对角线上排列 即构成特征值矩阵。用

12、eig函数计算的结果为：(2) p1 = -0.8483 -0.3290 0.41490.4514 -0.8589 0.24190.2767 0.3925 0.8771 d1 = -3.0645 0 00 1.7042 00 0 23.3603 两者的差别仅仅是特征值和特征向量的排列不同，因为eig 函数中把特征值按递增排列。 （3）det(A) = -122, det(d) = -122.0000, trace(A) = 22, trace(d) = 22, 可见其特征根的和仍为原矩阵的迹，其特征根的积仍为原 矩阵的行列式。 【例5-5-9】对称实矩阵对角化设有对称实矩阵试求 解：建模：矩阵

13、指数只有在A是对角矩阵时才可按对角 元素直接计算， 即： 注意它的非对角元素不符合指数运算规则， 。实矩阵对角化求矩阵指数(1)要利用这个性质，首先把A对角化。使A=pDp-1，将此式 代入y的表示式中，可得：为求其特征根和特征向量矩阵D, p。编成程序exn556： A = -2, 4, 1; 4, -2, -1; 1, -1, -3 , p,D = eig(A) 得p = -0.6572 -0.2610 0.70710.6572 0.2610 0.70710.3690 -0.9294 0.0000 D = -6.5616 0 00 -2.4384 00 0 2.0000实矩阵对角化求矩阵指

14、数(2) 由此得到： 算式中的最后一步可以用下列语句来实现。 y=p*exp(D(1,1),0,0;0,exp(D(2,2),0;0,0,exp(D(3,3)*inv( p) 得到 应用篇中与本节相关的实例(1) 【例6-1-1】 的物理数据的处理，用一根直线去拟合 有测量误差的数据点，这就是解一个超定的线性方程组 。每个方程都有误差，找到的最佳结果是使各方程的误 差均方值为最小，也就是最小二乘问题。 【例6-2-3】和【例7-1-2】的静力学平衡问题，因为 涉及两个物体，平衡方程和待求变量会超过四个。这些 通常都是线性方程组，用矩阵来解可以一次解出全部未 知数，最为方便快捷。 【例7-2-1

15、】的材料力学静不定问题和上述静力学平 衡问题相仿，只不过加了几个变形协调方程，这些方程 仍然是线性方程，所以不管方程和变量增加了多少，都 只要用一个矩阵方程一次就解出全部结果。只要系数输 入没有错，结果肯定是对的。应用篇中与本节相关的实例(2) 【例7-3-3】的二多自由度机械振动问题更是用矩阵解高阶微分 方程的典型。它不仅涉及代数方程的解，而且涉及特征值和特征向 量的计算问题。 【例8-1-1】的直流电路分析是典型的线性代数方程组，【例8-1 -5】的交流电路分析也化成线性代数方程组，只不过方程的系数矩 阵和未知数向量都是复数。在解复数方程组时，MATLAB更显示了 比手工解题的极大优越性，包括绘制相量图。【例8-1-8】网络参 数都是用矩阵表示和矩阵运算的。 【例8-2-2】的晶体管放大电路分析得出了一个复数矩阵方程， 它的所有系数元素又都是频率的函数。给了不同的频