4数值计算

1MATLAB基础与应用北京化工大学信息科学与技术学院主讲教师：郭青[e-mail]第4章数值计算2主要内容•4.1矩阵的基本运算•4.2矩阵的基本函数运算•4.3函数的零点•4.4函数的极值点•4.5数值积分•4.6常微分方程34.1矩阵的基本运算与数运算基本算术运算a=1.000020.800040.600060.400080.2000100.0000?a-1ans=019.800039.600059.400079.200099.0000?a*2ans=2.000041.600081.2000120.8000160.4000200.0000点积计算点乘，也叫向量的内积、数量积。指两个向量在其中一个向量方向上的投影的乘积。两个向量a=[a1,a2,…,an]和b=[b1,b2,…,bn]的点积定义为:函数dot(a,b)a,b必须同维。a=[123];b=[3,4,5];dot(a,b)ans=26sum(a.*b)%sum求和函数ans=26njjjbaba1)(向量的叉积叉积表示过两相交向量的交点的垂直于两向量所在平面的向量。cross（a，b）a，b必须为三维向量。混合积?c=cross(a,b)c=-24-2?dot(a,cross(b,c))ans=24矩阵的基本运算加减运算要求两矩阵必须同阶。?a=[123;234;345];?b=[111;222;333];?c=a+bc=234456678乘法要求a为i×j阶，b为j×k阶时，ab才能相乘。?e=[b,[555]']e=111522253335?f=a*ef=141414302020204526262660除法左除“\”:相当于Ax=B的解，x=A-1B。右除“/”：相当于xA=B的解，x=BA-1A-1B=(B’A’-1)’。通常，右除稍快一些，而左除可以避免奇异性。对于Ax＝B，其中A为（n×m）阶矩阵：n=m且非奇异时，方程为恰定方程；nm方程为超定方程；nm方程为欠定方程。矩阵的转置用标号’表示例：求矩阵的转置A=magic(4)%创建一个4阶魔方矩阵B=A’%求矩阵A的转置C=A(:,1)’%计算列向量的转置，得到行向量?A=[123;456;780;135];?B=[135;246];?A/Bans=00.5000-3.00003.5000-12.000010.25001.00000.0000?(B'\A')'ans=00.5000-3.00003.5000-12.000010.25001.00000.0000矩阵的逆运算•逆矩阵的定义：对于任意阶n×n方阵A，如果能找到一个同阶的方阵V，使得满足：A*V=I。其中I为n阶的单位矩阵eye(n)。则V就是A的逆矩阵。数学符号表示为：V=A-1。逆矩阵V存在的条件是A的行列式不等于0。•格式：V=inv(A)功能：返回方阵A的逆矩阵V?A=[21-3-1;3107;-124-2;10-15];?inv(A)ans=-0.04710.5882-0.2706-0.94120.3882-0.35290.48240.7647-0.22350.2941-0.0353-0.4706-0.0353-0.05880.04710.2941说明：如果矩阵不可逆，将给出错误信息也可以用表达式A^(-1)求得矩阵的行列式运算计算方阵A的行列式值函数det?A=[21-3-1;3107;-124-2;10-15];?a1=det(A)a1=-85?a2=det(inv(A))a2=-0.0118?a1*a2ans=1伪逆矩阵函数pinv伪逆矩阵的MATLAB定义：从数学意义上讲，当矩阵A为非方阵时，其矩阵的逆是不存在的。在MATLAB中，为了求线性方程组的需要，把inv(A′*A)*A′的运算定义为伪逆函数pinv，这样对非方阵，利用伪逆函数pinv可以求得矩阵的伪逆，伪逆在一定程度上代表着矩阵的逆。格式：C=pinv(A)功能：计算非方阵A的伪逆矩阵。矩阵的幂运算与数字的幂运算形式相同，用“^”算符。矩阵的指数运算常用函数expmexpm1expm2expm3矩阵的对数运算函数logm矩阵的开方运算函数sqrtm?b=magic(3)b=816357492?sqrtm(b)ans=2.7065+0.0601i0.0185+0.5347i1.1480-0.5948i0.4703+0.0829i2.0288+0.7378i1.3739-0.8207i0.6962-0.1430i1.8257-1.2725i1.3511+1.4155i?b^0.5ans=2.7065+0.0601i0.0185+0.5347i1.1480-0.5948i0.4703+0.0829i2.0288+0.7378i1.3739-0.8207i0.6962-0.1430i1.8257-1.2725i1.3511+1.4155i4.2矩阵的基本函数运算特征值与特征向量设n阶方阵A，如果数λ和n维非零列向量x使关系式Ax=λx成立，则λ称为A的特征值，x成为A对应于λ的特征向量函数[x,y]=eig(A)可以给出特征值和特征向量的值x为特征向量矩阵，y为特征值矩阵。?A=[73-2;34-1;-2-13];?[x,y]=eig(A)x=0.57740.0988-0.8105-0.5774-0.6525-0.49080.5774-0.75130.3197y=2.00000002.39440009.60564.2矩阵的基本函数运算矩阵的迹矩阵的迹指矩阵的对角线元素之和，也等于矩阵的特征值的和。函数traceA=[121;212;112];a=trace(A)a=4E=eig(A)E=4.3028-1.00000.6972b=sum(E)b=4.0000说明：E=eig(A)，求矩阵A的特征值列向量b=sum(E)，求矩阵E的列元素向量之和4.2矩阵的基本函数运算矩阵翻转函数fliplrflipudrot90a=73-234-1-2-13?fliplr(a)ans=-237-1433-1-2?flipud(a)ans=-2-1334-173-2?rot90(a)ans=-2-1334-173-2矩阵的秩e=111522223335?rank(e)ans=2迹函数矩阵所有对角线上元素的和称为矩阵的迹。函数trace正交空间函数函数orth用来求矩阵的一组正交基。条件数函数判断矩阵的“病态”程度。函数cond计算矩阵的条件数的值condest计算矩阵的1范数条件数的估计值rcond计算矩阵的条件数的倒数值伪逆函数函数pinv求解“病态”问题时，避免产生伪解。?a=magic(4)a=16231351110897612414151?b=a*[1111]';?inv(a)*bWarning:Matrixisclosetosingularorbadlyscaled.Resultsmaybeinaccurate.RCOND=1.567374e-017.ans=0800?pinv(a)*bans=1.00001.00001.00001.0000写成矩阵形式可表示为：AX＝B或XA＝B。其中系数矩阵A的阶数为m×n。一般线性方程组4.2.2线性代数方程求解mnmnmmnnnnbxaxaxabxaxaxabxaxaxa............22112222212111212111其中mnmmnnaaaaaaaaaA.........212222111211nxxxX21mbbbB21.,0)det(方程组存在唯一解约定ABAX线性方程组)det()det(AAxCramerii法则：引入矩阵除法求解。(1)求解方程AX=B格式：X=A\B条件：矩阵A与矩阵B的行数必须相等。(2)求解方程XA=B格式：X=B/A条件：矩阵A与矩阵B的列数必须相等。mnmnmmnnnnbxaxaxabxaxaxabxaxaxa12211212222121111212111一般线性方程组4.2.2线性代数方程求解3)如m=n则：X=A-1B即：X=inv(A)*B矩阵的分解三角(LU)分解函数lu所谓三角解就是将一个方阵表示成两个基本三角阵的乘积(A=LU），其中一个为下三角矩阵L，另一个为上三角形矩阵U，因而矩阵的三角分解又叫LU分解或叫LR分解。矩阵分解的两个矩阵分别可表示为：nnijaA}{100100012121nnlllLnnnnuuuuuuU000222112114.2.2线性代数方程求解LU分解主要调用格式：•[L，U]=lu(X)，其中Ｘ是任意方阵，L是下三角矩阵，U是上三角矩阵，满足的条件式为；•[L，U，P]=lu(X)，其中Ｘ是任意方阵，L是下三角矩阵，U是上三角矩阵，Ｐ是置换矩阵，满足的条件式为PX=LU。•Ｙ=lu(X)，其中Ｘ是任意方阵，把下三角矩阵和上三角矩阵合并在矩阵Ｙ中给出，满足的条件式为Y=L+U-I，该命令将损失置换矩阵Ｐ的信息。注意：进行lu分解时，矩阵X必须是方阵。4.2.2线性方程组求解【例】已知矩阵A=[5-12;12-1;314]，求其LU分解。A=[5-12;12-1;314];[L，U]=lu(A)L*U%测试分解的正确性4.2.3恰定方程组的解对于方程Ax=b，A为n×m矩阵，有三种情况：•当n=m时，此方程成为“恰定”方程；•当nm时，此方程成为“超定”方程；•当nm时，此方程成为“欠定”方程。4.2.3恰定方程组的解对于恰定方程组，最常用的方程组解法有：•x=inv(A)*b—采用求逆运算解方程；•x=A\b—采用左除运算解方程；•x=U\(L\b)—利用lu法求解。【例4-2】已知方程组，应用LU分解法求方程组的解。532413121215xA=[5-12;12-1;314];[L，U]=lu(A)b=[235]'x=U\(L\b)【例4-3】利用前２种方法对例4-2方程组进行求解。方法一：逆阵法A=[5-12;12-1;314];b=[235]';x=inv(A)*b方法二：左除法A=[5-12;12-1;314];b=[235]';x=A\b正交(QR)分解函数将矩阵A分解为一个正交矩阵与另一个矩阵的乘积称为矩阵A的正交分解。格式一：[Q,R]=qr(A)功能：产生与A同维的上三角矩阵R和一个实正交矩阵或复归一化矩阵Q，满足：A=Q*R，Q’*Q=I。格式二：[Q,R,E]=qr(A)功能：产生一个置换矩阵E，一个上三角矩阵R(其对角线元素降序排列)和一个归一化矩阵Q，满足A*E=Q*R；求解方法：X=R\(Q\B)4.3函数的零点•函数的零点是使f(x)=0的实数x，用二维图形表示也就是函数y=f(x)的图形和x轴交点的横坐标，即y=f(x)中y=0时x的取值。对于任意求解函数，可能有一个解，也可能有多个解，因此没有一个通用的解法。一般来说，可以通过作图，观察零点所在的位置，然后通过计算不断地逼近真实值，并在一定精度时停止计算。4.3.1一元函数的零点在MATLAB中求解一元函数的零点的函数是fzero()，其调用格式如下：•x=fzero(fun，x0)，参数fun为一元函数名，x0表示求解的初始数值。•[x,fval,exitflag,output]=fzero(fun,x0,options)，options是优化迭代所采用的参数选项，在输出参数中，fval表示对应的函数值，exitflag表示程序退出的类型，output为优化