1、线性代数在数学建模中的应用第六章 线性代数模型 6.1 Matlab求解线性代数工具简介矩阵的秩. rref或rrefmovie格式R=rref(A)%用高斯约当消元法和行主元法求A 的行最简行矩阵R.rrefmovie(A)%给出每一步化简的过程.方阵的行列式:det(A)逆矩阵:inv(A)指令inv(A)给出方阵A的逆矩阵,如果A不可逆,则inv(A)给出的矩阵的元素都是Inf 利用初等变换也可以求出逆矩阵,构造n行2n列的矩阵(A E),并进行行初等变换,当把A变为单位矩阵时,E就变成了A的逆矩阵利用matlab命令rref可以求出矩阵的行简化阶梯形输入命令:D=A,eye(3) D
2、= 1 2 3 1 0 0 2 2 1 0 1 0 3 4 3 0 0 1 rref(D) ans = 1.0000 0 0 1.0000 3.0000 -2.0000 0 1.0000 0 -1.5000 -3.0000 2.5000 0 0 1.0000 1.0000 1.0000 -1.0000 线性方程组的求解是用矩阵除来完成的,,当且可逆时,给出唯一解这时矩阵除相当于;当时,矩阵除给出方程的最小二乘解;当时,矩阵除给出方程的最小范数解例6.1:解方程组:解:输入命令: a=1 -1 1 2;1 1 -2 1;1 1 1 0;1 0 1 -1; b=1;1;2;1; x=ab x =
3、0.8333 0.7500 0.4167 0.2500 输入命令: z=inv(a)*b z = 0.8333 0.7500 0.41670.2500 例6.2:解方程组:解:方程的个数和未知数不相等,用消去法,将增广矩阵化为行简化阶梯形,如果系数矩阵的秩不等于增广矩阵的秩,则方程组无解;如果系数矩阵的秩等于增广矩阵的秩,则方程组有解,方程组的解就是行简化阶梯形所对应的方程组的解输入命令: a=2 1 1 -1 -2 2;1 -1 2 1 -1 4;2 -3 4 3 -1 8; rref(a) ans = 1 0 0 0 0 0 0 1 0 -1 -1 0 0 0 1 0 -1 2 由结果看出
4、,,为自由未知量,方程组的解为: 例:解方程组:解:输入命令: a=1 -1 -1 1;1 -1 1 -3;1 -1 0 -1;1 -1 -2 3; rref(a) ans = 1 -1 0 -1 0 0 1 -2 0 0 0 0 0 0 0 0 由结果看出,,为自由未知量,方程组的解为: 6.3 交通流量模型城市道路网中每条道路、每个交叉路口的车流量调查,是分析、评价及改善城市交通状况的基础。根据实际车流量信息可以设计流量控制方案,必要时设置单行线,以免大量车辆长时间拥堵。下面考虑一简化的交通流量统计问题。问题描述:某城市有下图所示的交通图,每条道路都是单行线,需要调查每条道路每小时的车流量
5、。图中的数字表示该条路段的车流数。如果每个交叉路口进入和离开的车数相等,整个图中进入和离开的车数相等。试解决以下问题:(1)建立确定每条道路流量的线性方程组; (2)分析哪些流量数据是多余的;(3)为了唯一确定未知流量, 需要增添哪几条道路的流量统计。图6.1 某城市单行线车流量模型假设:根据题目要求,作以下假设:(1)全部流入网络的流量等于全部流出网络的流量,这个从图中信息可以确定;(2)全部流入一个节点的流量等于全部流出此节点的流量。模型建立: 由网络流量假设,每条道路流量的线性方程组为即 (6.3)模型求解:上述方程组的系数矩阵为在Matlab中利用命令rref(B)将增广矩阵化为阶梯形
6、行最简形式其对应的非齐次方程组可表示为 (6.4)取(x9, x10, x11, x12)为自由未知量,分别赋四组值为(1,0,0,0),(0,1,0,0),(0,0,1,0),(0,0,0,1),则(6.4)对应齐次方程组基础解系中四个解向量为 取自由未知量(x9, x10, x11, x12)为(0,0,0,0)得非齐次方程组(6.4)的特解于是方程组(6.3)的通解为其中k1,k2,k3,k4为任意常数,x的每一个分量即为交通网络未知部分的具体流量,显然它有无穷多解。结果分析:从增广矩阵的阶梯形行最简形式可以看出,方程组(6.3)中最后一个方程是多余的,即最后一个方程中的流量数据150和
7、290不用统计。而如果要唯一确定未知流量,可以增添x9, x10, x11, x12(或x2, x9, x10, x12等)这几条道路的流量统计。 在特解中,x6,x7为负,这表明对应的两条单行线应该改变方向才合理。6.4 化学方程式配平问题在用化学方法处理污水过程中, 有时会涉及到复杂的化学反应。这些反应的化学方程式是分析计算和工艺设计的重要依据。在定性地检测出反应物和生成物之后,可以通过求解线性方程组配平化学方程式。问题描述:某厂废水中含KCN, 其浓度为650mg/L。现用氯氧化法处理, 发生如下反应: KCN + 2KOH + Cl2 = KOCN + 2KCl + H2O.投入过量液
8、氯, 可将氰酸盐进一步氧化为氮气。请配平下列化学方程式: KOCN + KOH + Cl2 = CO2 + N2 + KCl + H2O. (注: 题目摘自福建省厦门外国语学校2008-2009学年高三第三次月考化学试卷)模型建立:设x1KOCN x2KOH x3Cl2 = x4CO2 x5N2 x6KCl x7H2O,则根据平衡关系有即模型求解:在Matlab中利用命令rref(B)将增广矩阵化为阶梯形行最简形式A = 1,1,0,0,0,-1,0; 1,1,0,-2,0,0,-1; 1,0,0,-1,0,0,0; 1,0,0,0,-2,0,0; 0,1,0,0,0,0,-2; 0,0,2,
9、0,0,-1,0;x = null(A,r); %求齐次方程组的基础解系format rat, x %以有理数形式输出Matlab执行后得ans = 1 2 3/2 1 1/2 3 1可见上述齐次线性方程组的通解为x = k(1, 2, 3/2, 1, 1/2, 3, 1)T。 取k = 2得x = (2, 4, 3, 2, 1, 6, 2)T。 可见配平后的化学方程式如下2KOCN + 4KOH + 3Cl2 = 2CO2 + N2 + 6KCl + 2H2O. 小结:利用线性方程组配平化学方程式是一种待定系数法。关键是根据化学方程式两边所涉及到的各种元素的量相等的原则列出方程。所得到的齐次
10、线性方程组Ax =b中所含方程的个数等于化学方程式中元素的种数s, 未知数的个数就是化学方程式中的项数n。 当r(A) = n 1时, Ax =b的基础解系中含有1个(线性无关的)解向量。这时在通解中取常数k为各分量分母的最小公倍数即可。 例如本例中1, 2, 3/2, 1, 1/2, 3, 1分母的最小公倍数为2, 故取k = 2。当r(A) n 2时, Ax =b的基础解系中含有2个以上的线性无关的解向量。这时可以根据化学方程式中元素的化合价的上升与下降的情况, 在原线性方程组中添加新的方程。6.6小行星的轨道模型问题描述: 一天文学家要确定一颗小行星绕太阳运行的轨道,他在轨道平面内建立以
11、太阳为原点的直角坐标系,在两坐标轴上取天文测量单位(一天文单位为地球到太阳的平均距离:9300万里)。在5个不同的时间对小行星作了5次观察,测得轨道上5个点的坐标数据如表6.1.表6.1 坐标数据x1x2x3x4x5X坐标4.55965.08165.55465.96366.2756y1y2y3y4y5Y坐标0.81451.36851.98952.69253.5265由Kepler(开普勒)第一定律知,小行星轨道为一椭圆。试确定椭圆的方程以及各个参数,以供研究。问题分析与模型建立: 由Kepler第一定律知,小行星轨道为一椭圆。而椭圆属于二次曲线,二次曲线的一般方程为. (6.6)天文学家确定小
12、行星运动的轨道,依据的是轨道上五个点的坐标数据:(x1, y1),(x2, y2),(x3, y3),(x4, y4),(x5, y5)。为了确定方程中的五个待定系数,将五个点的坐标分别代入上面的方程,得这是一个包含五个未知数的线性方程组,写成矩阵利用题目数据,可以求解这一线性方程组,进而确定方程(6.6),但所得的是一个二次曲线方程的一般形式。为了知道小行星轨道的一些参数,还必须将二次曲线方程化为椭圆的标准方程形式:由于太阳的位置是小行星轨道的一个焦点,这时可以根据椭圆的长半轴和短半轴计算出小行星的近日点和远日点距离。根据解析几何中二次曲线理论,椭圆经过旋转和平移两种变换后方程可表示为:,这
13、里, 其中D为二次曲线的矩阵,为矩阵C的特征根(不妨设),即二次曲线(6.6)的特征根。所以,椭圆长半轴:;椭圆短半轴:;椭圆半焦距:。而椭圆的中心坐标满足下面的方程组:因为原点即太阳是椭圆的一个焦点,因此利用两个焦点关于中心对称,可知另一个焦点的坐标为。模型求解:根据以上分析,在Matlab中编写如下程序:format shortclear;clc;x=4.5596 5.0816 5.5546 5.9636 6.2756;y=0.8145 1.3685 1.9895 2.6925 3.5265;plot(x,y,r*)hold on%求解线性方程组,确定二次曲线的系数b=-ones(leng
14、th(x),1);A=x.2 2*x.*y y.2 2*x 2*y;a=Ab %绘制小行星轨道椭圆曲线syms s t eq1=a(1)*s2+2*a(2)*s*t+a(3)*t2+2*a(4)*s+2*a(5)*t+1;ezplot(eq1,-2,8,-2,8) title(小行星轨迹) hold ongrid on%小行星轨道的相关参数C=a(1) a(2);a(2) a(3);D=a(1) a(2) a(4);a(2) a(3) a(5);a(4) a(5) 1;X=C-a(4);-a(5); %确定椭圆中心坐标x0=X(1);y0=X(2);U,d=eig(C); %求C的特征值和特征
15、向量 f=det(D)/det(C); %计算C和D的行列式的比值a0=sqrt(-f/d(2,2) %确定椭圆半长轴和半短轴b0=sqrt(-f/d(1,1)c=sqrt(a02-b02) %确定椭圆半焦距h=a0-c %确定小行星近日点距离H=a0+c %确定小行星远日点距离x1=2*x0;y1=2*y0; %确定椭圆除原点外另一焦点坐标plot(x0,y0,ro,0,0,b+,x1,y1,b+) %在图中标注中心和两个焦点%动画演示小行星运动轨迹t=0:0.1:2*pi;u=a0*cos(t);v=b0*sin(t);U=U(:,2),U(:,1);V=U*u;v;x1=V(1,:)+x0;y1=V(2,:)+y0;n=length(t);s=1:n;si=1:0.1:n;x2=interp1(s,x1,si);x2=x2,x2,x2,x2,x2;y2=interp1(s,y1,si);y2=y2,y2,y2,y2,y2;comet(x2,y2)运行后整理得到以下参数值:二次曲线的系数: 半长轴a0 =4.3799;半短轴b0 =2.4299;焦距c =3.6440;近日点距离h =0.7359;远日点距离H=8.0239.椭圆轨迹见图6.1图6.1 小行星运行椭圆轨迹图中*表示测量数据点,o代表椭圆中心,+代表两个焦点。
copyright@ 2008-2022 冰豆网网站版权所有
经营许可证编号:鄂ICP备2022015515号-1