基于MATLAB的数值分析

基于MATLAB的数值分析以软件MATLAB作为辅助工具介绍数值分析（科学与工程计算）的基本内容，注重讲授一些求解方程以及结果可视化的知识和技巧，使同学们能够有效地解决问题并处理计算结果。内容包括：1、MATLAB编程和绘图2、数值分析的数学基础3、数值算法在工程、科学中的应用第一章MATLAB入门1、MATLAB的命令窗口工作都在此处完成。2、怎样进行计算运算对象：矩阵算术运算符：+-*^.*./.^倒除：\a\b=b/a变量与变量名：变量名和变量名类型不需声明。3、数据显示格式默认格式：5位（formatshort)formatlong16位formate短的浮点格式formatlonge长的浮点格式4、清除命令clear:清除所有使用过的变量或某个（些）变量clc:清除命令窗口5、程序结构分支：if__else__end;if__elseif__end;if__break__end循环：for__end;while__end6、读写输入数据：z=input(‘typeyoueinput:’)键盘输入格式化输出：fprintf(‘e_format%12.5e\n’,vol)7、数学函数8、功能函数sort(x)sum(x)max(x)min(x)mod(x,y)rand(n)eval(s)9、编程（编写M文件）10、绘图第二章数值代数内容：数值代数就是研究有关矩阵计算的问题。主要包括：1、线性代数方程组的求解；2、矩阵特征值问题要求：1、掌握用MATLAB求解的方法2、知道那些问题是困难的，那些问题是不可解的。A=zeros(m,n)m行n列的零矩阵I=eye(n)n阶单位矩阵A=ones(m,n)元素均为1A’A的转置A(:,k)A(k,:)A(m1:m2,n1:n2)inv(A)A的逆size(A)A的大小hilb(n)Hilbert矩阵2.1矩阵bAx情况1：m=n(正规方程），最常见；情况2：mn(不定方程）；情况3：mn(超定方程）；本节只介绍情况1。倍计算时间大约是上面的效率最高50*)(\bAinvxbAxMATLAB命令：2.2解线性代数方程组的MATLAB命令线性代数方程组并不总是数值可解的。只有当矩阵A的行列式不为零时才行！矩阵A的行列式即使不为零，但当很小或很大时，解的误差可能很大。计算矩阵行列式的MATLAB命令：)det(AD2.4病态问题有许多线性代数方程组理论上是可解的，但实际计算中由于受到舍入误差的影响而无法得到精确解。此类问题成为病态问题。病态问题的计算过程中，小的舍入误差或系数矩阵的微小变化都可能使解产生很大误差。（例子P97)2.3不可解问题病态矩阵的一个重要标志是条件数：1)(AAAcondMATLAB命令：)(Acond当矩阵是病态时，其条件数一定很大，但它并不能直接说明解的误差。线性方程组解的误差程度也取决于计算环境的精度。条件数和行列式与计算环境是相互独立的。所以大条件数或小行列式未必意味无法直接精确求得线性方程组的解，它只意味着有很大误差可能。而实际上如果采用更高精度的计算环境则很可能得到非常满意的解。Hilbert矩阵是非常著名的病态矩阵（hilb(n)),它经常用来检验算法的数值稳定性的好坏。两种原因使我们想了解求解线性代数方程组的算法。一是实际工作中要用其它计算机语言（Fortran&C等）编写应用程序；二是MATLAB处理大型稀疏矩阵方程组显得很笨拙或无能为力。nnnnnbabaabaaabABbAx2222111211()(初等行变换）2.5线性代数方程组的求解方法（算法）由线性代数的理论：下面讨论实现过程：第一步：消元。进行到第k步时必有nnnknnkknkkkkkkknkkkknkkbaaabaaabaaabaaaa1,1,11,1,11,111,1111a(k,k)作为主元，第k行依次乘a(i,k)/a(k,k)加到第i行,i=k+1:n。总共n-1步完成，k=1:n-1.一、Gauss消去法endendendendjkakiajiajia:nkjforkbkiaibibkkakiakiankiforkkaifnkfor);,(),(),(),(1);(),()()();,(/),(),(:1end;break,,0),(1:1当a(k,k)=0，则上述消去法无法进行；或当其绝对值相对太小可能会出现大的计算误差。选主元法可避免这种情况。下面介绍常用的按列选主元的Gauss法。endendendjkakiajiajia:nkjforkbkiaibibkkakiakiankifortbkkbtankkkakkbkbnkkkankkab(k)tba(k,k:n);takkkkkknkaabskkqnkfor);,(),(),(),(1);(),()()();,(/),(),(:1;)(;):,();()();:,():,(;;1)));,:((max(],[1:1列主元Gauss法运算量(只考虑乘除运算）：第k步=n-k次除法+n-k次乘法+（n-k)*(n-k)次乘法;)(31)(2)(32311112nonnnknknnknk总的乘除运算量=第二步：回代求解1:1);,(/)),(();,(/1niiiaxjiabxnnabxnijjiinnendiiadibibendjbjiaddnijfordnifornnanbnb);,(/))(()();(),(:1;01:1:1);,(/)()(%%%%%%%%二、LU分解法LU分解的目的是将矩阵A转换为两个矩阵的乘积，即。方程为上三角矩阵。此时，是单位下三角矩阵，yUxbLybLUxbAxLLUA,U,好处是：对于线性方程组，如果需要多次求解不同的非齐次项，此时LU分解的效率将大大超过高斯消去法。LU分解的MATLAB命令：[l,u]=lu(A)和[l,u,p]=lu(A)前面讲到的不选主元的高斯消去法和列主元高斯消去法将能实现LU分解。不选主元的高斯消去法用于下面两类矩阵肯定能成，即严格对角占优矩阵或对称正定矩阵，其他矩阵就难说了！列主元高斯消去法是解决一般中小型稠密矩阵方程组最有效的方法之一。下面讲解列主元高斯消去法实现LU分解的算法。.1:1))).,:((max(],[,1)(.1),(),(1),(),1(11i)jP(i,nkknkaabskkqkkkkkkkkkkaknakkakkaLjk为其中行交换后的矩阵，行与第为单位矩阵第记1、LU分解的代数理论)).(,(,333223221131211112222113333223221131211112223322221121111kkkkPPUaaaaaaaaaaAPLPLPLPLaaaaaaaaaaaAPLPLaaaaaaaaaAPLknnnnnnnnnnnnnnnnnnnnn其中，于是：..~,,~,,~,~;)()()())()((1211212111211211212111211212112112112321121111222211PPPPPPLPPLPPPLPPPLPLPLLLUAPPPPPLPPPPPLPPPPPLPPPLPLUAPLPLPLPLnnnnkkkkknknnnnnnnnnnnkkkkknnnnnnnnnnnnnn令：.~~~,.~~~111211121为单位下三角矩阵。可以证明其中LLLLLLUPAUPALLLnnnn现在我们只要将列主元高斯消去法稍加改造即是LU分解的算法。列主元高斯消去法的矩阵表示：.),,,),:1(U))(,()(.',2,111222211knkkkkkknnnnLaaaknkAAAkkkkPPkkkUAPLPLPLPL唯一确定了（部分，即的下三角的每一输出的矩阵；的上三角为输出的矩阵；唯一确定了2、LU分解算法endendendjkakiajiajia:nkjforkkakiakiankifortankkkkankkkkankkaa(k,k:n);tamkkkkknkaabsmqnkforasizenmaLGaussakkfunction);,(),(),(),(1);,(/),(),(:1;):),(();:),(():,(;1)()));,:((max(],[1:1);(],[)(],[运算量和高斯消去法一样endendt;,k)a(kk(i);a(kk(i),k)a(i,k)a(i,k);t:n-kifor:n-forkendt;k(k),:)p(k;p(kk(k),:)p(k,:)kptnkforneyepaLGaussakkasizenmajiesanjiaofenpulfunction1121:);,(1:1);();(],[);(],[)(],,[endenda(i,j);u(i,j)i:njfor:nforiendendenda(i,j);i,j)l(else;i,j)l(jifi:ijfor:nfori11113、三对角矩阵方程的追赶法三对角且主对角严格占优矩阵方程是一类来源丰富的问题。比如，微分方程数值解或样条插值等问题中的正规方程组。解这种问题必须考虑其矩阵稀疏的特征，减少算法的计算量。nnnabcabcaT12211三对角矩阵形如：fTx问题：求解T的LU分解具有形式：nnnucucuUllL12112,111由T=LU推得：).:2(,/;1111niclauublauiiiiiii1:1;/)(,/:2;,1111niuxcyxuyxniylfyfyyUxfLyfTxiiiiinnniiiiend1))/a(i);f(i*c(i)-(f(i)f(i)1:-1:1-niforf(n)/a(n);f(n)end1);-f(i*b(i)-f(i)f(i)n:2iforend1);-c(i*b(i)-a(i)a(i)1);-b(i)/a(ib(i)n:2iforlength(a);nf)c,b,a,zhuiganfa(ffunction最终解为f.乘除运算量：5n=O(n).4、对称正定矩阵方程的cholesky分解法对称正定矩阵方程的来源比较丰富，比如线性回归、拟合等问题。解这类问题必须考虑其矩阵对称的特征，减少算法的计算量。为对称正定矩阵。其中问题：求解A,bAx由于A为对称正定矩阵，A必有cholesky分解：nnnTdddlllLLDLA212121D,111,其中jjkjkikkijijjkkjkjjjiiTdlldaldladnidaladLDLA/)(n:1jin,:2j:2,/,11112111111得到关系：的两边，比较endendj);t)/a(j,-j)(a(i,j)a(i,endk);a(j,*k)a(i,*k)a(k,tt1-j:1kfor0;tn:1jifort;-