基于MATLAB的数值分析.ppt

基于MATLAB的数值分析,以软件MATLAB作为辅助工具介绍数值分析（科学与工程计算）的基本内容，注重讲授一些求解方程以及结果可视化的知识和技巧，使同学们能够有效地解决问题并处理计算结果。,内容包括：1、MATLAB编程和绘图2、数值分析的数学基础3、数值算法在工程、科学中的应用,第一章 MATLAB入门,1、MATLAB的命令窗口 工作都在此处完成。2、怎样进行计算 运算对象：矩阵 算术运算符：+-*.*./.倒除：ab=b/a 变量与变量名：变量名和变量名类型不需声明。,3、数据显示格式 默认格式：5位（format short)format long 16位 format e 短的浮点格式 format long e 长的浮点格式4、清除命令 clear:清除所有使用过的变量或某个（些）变量 clc:清除命令窗口5、程序结构 分支：if _else_end;if_elseif_end;if_break_end 循环：for_end;while_end,6、读写 输入数据：z=input(type youe input:)键盘输入 格式化输出：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)元素均为1 A A的转置 A(:,k)A(k,:)A(m1:m2,n1:n2)inv(A)A的逆 size(A)A的大小 hilb(n)Hilbert矩阵,2.1 矩阵,情况1：m=n(正规方程），最常见；情况2：mn(超定方程）；本节只介绍情况1。,MATLAB命令：,2.2 解线性代数方程组的MATLAB命令,线性代数方程组并不总是数值可解的。只有当矩阵A的行列式不为零时才行！矩阵A的行列式即使不为零，但当很小或很大时，解的误差可能很大。,计算矩阵行列式的MATLAB命令：,2.4 病态问题 有许多线性代数方程组理论上是可解的，但实际计算中由于受到舍入误差的影响而无法得到精确解。此类问题成为病态问题。病态问题的计算过程中，小的舍入误差或系数矩阵的微小变化都可能使解产生很大误差。（例子 P97),2.3 不可解问题,病态矩阵的一个重要标志是条件数：,MATLAB命令：,当矩阵是病态时，其条件数一定很大，但它并不能直接说明解的误差。线性方程组解的误差程度也取决于计算环境的精度。条件数和行列式与计算环境是相互独立的。所以大条件数或小行列式未必意味无法直接精确求得线性方程组的解，它只意味着有很大误差可能。而实际上如果采用更高精度的计算环境则很可能得到非常满意的解。Hilbert矩阵是非常著名的病态矩阵（hilb(n),它经常用来检验算法的数值稳定性的好坏。,两种原因使我们想了解求解线性代数方程组的算法。一是实际工作中要用其它计算机语言（Fortran&C等）编写应用程序；二是MATLAB处理大型稀疏矩阵方程组显得很笨拙或无能为力。,2.5 线性代数方程组的求解方法（算法）,由线性代数的理论：,下面讨论实现过程：第一步：消元。进行到第k步时必有,a(k,k)作为主元，第k行依次乘a(i,k)/a(k,k)加到第i行,i=k+1:n。总共n-1步完成，k=1:n-1.,一、Gauss消去法,当a(k,k)=0，则上述消去法无法进行；或当其绝对值相对太小可能会出现大的计算误差。选主元法可避免这种情况。下面介绍常用的按列选主元的Gauss法。,列主元Gauss法,运算量(只考虑乘除运算）：第k步=n-k次除法+n-k次乘法+（n-k)*(n-k)次乘法;,总的乘除运算量=,第二步：回代求解,%,二、LU分解法,LU分解的目的是将矩阵A转换为两个矩阵的乘积，即,好处是：对于线性方程组，如果需要多次求解不同的非齐次项，此时LU分解的效率将大大超过高斯消去法。,LU分解的MATLAB命令：l,u=lu(A)和l,u,p=lu(A),前面讲到的不选主元的高斯消去法和列主元高斯消去法将能实现LU分解。不选主元的高斯消去法用于下面两类矩阵肯定能成，即严格对角占优矩阵或对称正定矩阵，其他矩阵就难说了！列主元高斯消去法是解决一般中小型稠密矩阵方程组最有效的方法之一。下面讲解列主元高斯消去法实现LU分解的算法。,1、LU分解的代数理论,现在我们只要将列主元高斯消去法稍加改造即是LU分解的算法。列主元高斯消去法的矩阵表示：,2、LU分解算法,运算量和高斯消去法一样,3、三对角矩阵方程的追赶法,三对角且主对角严格占优矩阵方程是一类来源丰富的问题。比如，微分方程数值解或样条插值等问题中的正规方程组。解这种问题必须考虑其矩阵稀疏的特征，减少算法的计算量。,三对角矩阵形如：,T的LU分解具有形式：,由T=LU推得：,最终解为f.乘除运算量：5n=O(n).,4、对称正定矩阵方程的cholesky分解法,对称正定矩阵方程的来源比较丰富，比如线性回归、拟合等问题。解这类问题必须考虑其矩阵对称的特征，减少算法的计算量。,由于A为对称正定矩阵，A必有cholesky分解：,三、线性代数方程组的迭代法,线性代数方程组的迭代法并不适用于所有问题，但它对一些特定类型的问题非常有效。当问题是大型稀疏矩阵方程时，高斯消去法的效率会变得非常低，而且有时还会超出内存要求。对于这样的问题就需要使用迭代法。大系统问题的求解最终归结为大型稀疏矩阵方程。比如，电网络、场方程的数值计算、运筹问题等。尽管迭代法的种类很多，这里只介绍其中的三种：1、Jacobi迭代；2、Gauss-Seidel迭代;3、超松弛迭代法（SOR)。,1、线性代数方程组的迭代法的一般理论,迭代算法简单，但问题是：1、能否保证算法收敛？2、能否充分利用矩阵的稀疏性，使运算量和存贮量尽量少。,定理一 迭代法收敛的充分必要条件是,定理二 迭代法收敛的充分必要条件是,定理三 迭代法收敛的充分条件是,收敛性定理,对于任意的初值,迭代矩阵B的构造,充分利用矩阵的稀疏性，使运算量和存贮量尽量少的办法就是要求迭代矩阵B与原矩阵A有相同的稀疏结构。具体就是：,常用迭代法及其收敛性,定理4：SOR收敛的必要条件是02.定理5：如果A是严格对角占优矩阵或不可分弱对角占优矩阵，则1、Jacobi收敛；2、gauss-seidel收敛；3、当0=1 时，SOR必收敛。定理6：如果A是对称正定矩阵，则1、当2D-A正定时，Jacobi收敛；2、gauss-seidel收敛；3、当02时，SOR必收敛。其他矩阵不能保证收敛,四、线性代数方程组共轭梯度法,即，解对称正定矩阵方程组的问题等价于求二次函数极小点的问题。如何寻找f(x)的极小点呢？基本思想就是下降法。,1、下降法的基本原理,2、最速下降法,我们知道，函数的负梯度方向是函数值下降最快的方向。,function x=zuisu(a,b,x,eyta)n=length(b);r=b-a*x;while max(abs(r)eyta t=r*r/(r*a*r);x=x+t*r;r=b-a*x;end,这一算法有简单易行，可充分利用A的稀疏性等特点，但当A的最大特征值远远大于最小特征值时收敛速度变得非常慢，以至于完全不适用。最速下降法并非最速！下面的共轭梯度法可有效解决这一问题。,3、共轭梯度法,能否选出比负梯度方向更好的下降方向吗？能！,这一算法称为共轭梯度法，简称CG(conjugate gradient),实际计算中由于误差的影响，,之间的正交性很快就会,消失，以至于有限步内不能得到精确解，所以CG实际是一种迭代算法。,CG是保稀疏且便于并行处理的算法，它是求解对称正定矩阵方程最优秀的算法之一。而且目前它已推广到一般矩阵方程，称为预优共轭梯度算法。,实用共轭梯度算法如下：,function x=cg(a,b,x,eyta)n=length(b);r=b-a*x;p=r;q0=r*r;while q0eyta w=a*p;t=q0/(p*w);x=x+t*p;r=r-t*w;q=r*r;s=q/q0;p=r+s*p;q0=q;end,运算量：每一步迭代的乘除运算量为,2.6 矩阵特征值问题,一、幂法,MATLAB命令：D=eig(A)和X,D=eig(A),实用算法：,function lmta,x=mifa(a,epsl)m,n=size(a);x0=ones(n,1);while 1 x=a*x0;lmta,m=max(abs(x);lmta=sign(x(m)*lmta;x=x/lmta;if abs(x-x0)epsl,break,end;x0=x;end,一、幂法,二、反幂法及其原点位移,反幂法用来求A的按模最小的特征值。,实用算法：,function lmta,x=fanmifa(a,epsl)m,n=size(a);x0=ones(n,1);while 1 x=ax0;lmta,m=max(abs(x);lmta=sign(x(m)*lmta;x=x/lmta;if max(abs(x-x0)epsl,break,end;x0=x;endlmta=1/lmta;,1、反幂法,2、带原点位移的反幂法,反幂法与“原点位移”相配合，求指定点附近的某个特征值和特征向量，并可用于加速幂法的收敛性。,另一方面，s离A的某个特征值越近，收敛越快。因此不论用幂法求A的按模最大特征值，还是利用反幂法求A的按模最小特征值，为了加快收敛，均可以用迭代m步后的近似值lmta作为最初始的位移值,实行动态位移迭代。,动态位移幂法,function lmta,x=dongtamifa(a,epsl0,epsl)lmta,x=mifa(a,epsl0);x0=x;lmta0=lmta;n=length(x0);while 1 x=(a-lmta0*eye(n)x0;lmta,m=max(abs(x);lmta=sign(x(m)*lmta;x=x/lmta;lmta=1/lmta+lmta0;if max(abs(x-x0)epsl,break,end;x0=x;lmta0=lmtaend,动态位移反幂法,function lmta,x=dongtaifanmifa(a,epsl0,epsl)lmta,x=fanmifa(a,epsl0);x0=x;lmta0=lmta;n=length(x0);while 1 x=(a-lmta0*eye(n)x0;lmta,m=max(abs(x);lmta=sign(x(m)*lmta;x=x/lmta;lmta=1/lmta+lmta0;if max(abs(x-x0)epsl,break,end;x0=x;lmta0=lmta;end,三、幂法综述,1、幂法和反幂法只能用于求解可对角化矩阵的实数特征值和特征向量，不能求解复数特征值；2、幂法和反幂法均为线性收敛，收敛速度由收敛因子决定，效率不高。3、动态位移可以大大减小收敛因子加速收敛；4、不适于求解全部特征值。5、对称矩阵自然适用于幂法，此时采用2-范数标准化向量的算法至少平方收敛。,求解一般矩阵全部特征值的办法是下面的相似变换法,四、矩阵全部特征值的QR迭代算法,求解一般矩阵全部特征值的办法是相似变换法。其理论根据是：,1、矩阵的两种正交变换,平面旋转变换,镜面反射矩阵,镜面反射矩阵的意义是“成批”消去向量的非零元素。,function a,b,x=householder(x)%x gived,seek a househouder convert H=I-bxx,make Hx=-a.n=length(x);t=max(abs(x);if t=0 a=0;b=0;else x=x/t;a=sqrt(x*x);if x(1)0,a=-a;end x(1)=x(1)+a;b=1/(a*x(1);a=t*a;end,2、方阵的正交三角化分解,方阵的正交三角化分解,即,function arfa,bata,a=qrfenjie(a)m,n=size(a);for k=1:n-1 arfa(k),bata(k),a(k:n,k)=householder(a(k:n,k);for j=k+1:n a(k:n,j)=a(k:n,j)-bata(k)*a(k:n,k)*a(k:n,j)*a(k:n,k);endend,function q,r=qrgouzao(a)%实现a=qr,q是正交矩阵,r是上三角矩阵。m,n=size(a);arfa,bata,a=qrfenjie(a);q=eye(n)-bata(1)*a(1:n,1)*a(1:n,1);for k=2:n-1 q(1:k-1,k:n)=q(1:k-1,k:n)-bata(k)*q(1:k-1,k:n)*a(k:n,k)*a(k:n,k);q(k:n,k:n)=q(k:n,k:n)-bata(k)*q(k:n,k:n)*a(k:n,k)*a(k:n,k);endfor i=1:n-1 for j=i+1:n r(i,j)=a(i,j);endendfor i=1:n-1 r(i,i)=-arfa(i);endr(n,n)=a(n,n);,3、化矩阵为Hessenberg型,function q,h=hessenberghua(a);%h=qaq,h为Hessenberg矩阵，q为正交矩阵。m,n=size(a);q=eye(n);for k=1:n-2 arfa(k),bata(k),a(k+1:n,k)=householder(a(k+1:n,k);p=eye(n-k)-bata(k)*a(k+1:n,k)*a(k+1:n,k);a(1:k,k+1:n)=a(1:k,k+1:n)*p;a(k+1:n,k+1:n)=p*a(k+1:n,k+1:n)*p;q(1:n,k+1:n)=q(1:n,k+1:n)*p;endh=a;for i=1:n-2 h(i+1,i)=-arfa(i);endfor j=1:n-2 for i=j+2:n h(i,j)=0;endend,注：当A为对称矩阵时，其Hessenberg 形为对称三对角矩阵。,4、Hessenberg矩阵在QR分解下的不变性,显然，实现Hessenberg矩阵的QR分解用Givens变换合适。即执行n-1次Givens变换:,function q,h=givensqr(h)%h为不可约hessenberg矩阵，用givens变换进行QR分解。m,n=size(h);q=eye(n);for k=1:n-1 if h(k+1,k)=0 c=h(k,k)/sqrt(h(k,k)2+h(k+1,k)2);s=h(k+1,k)/sqrt(h(k,k)2+h(k+1,k)2);d=c,s;-s,c;h(k:k+1,k:n)=d*h(k:k+1,k:n);if k=1 q=d zeros(2,n-2);zeros(2,n-2)eye(n-2);else q=q*eye(k-1)zeros(k-1,2)zeros(k-1,n-k-1).;zeros(k-1,2)d zeros(2,n-k-1).;zeros(k-1,n-k-1)zeros(2,n-k-1)eye(n-k-1);end end end,5、基本QR算法,一般QR 迭代算法的收敛性比较复杂，这里不再介绍。仅指出，在一定条件下由基本QR算法生成的系列,收敛为准上三角矩阵，对角块按特征值的模从大,到小排列。,function q,h=jibenQR(a,epsl)m,n=size(a);q0,h=hessenberghua(a);t=zeros(1,n-1);while 1 for i=1:n-1 if abs(h(i+1,i)=(abs(h(i,i)+abs(h(i+1,i+1)*epsl h(i+1,i)=0;t(i)=i+1;end end k=1;while k=n-2 if t(k)=t(k+1)yes=1;k=n-1;else k=k+1;yes=0;end end if yes=0,break,end;q,h=givensqr(h);h=h*q;q=q0*q;q0=q;end,6、QR 算法的改善,基本QR算法计算量和存储量都很大，且若收敛是线性的。因此，需要从这两方面加以改善。,（1）划分和收缩（2）原点位移,第三章 数值逼近,内容：数值逼近也称数据模拟，是为离散数据（不论何种途径得到）建立简单连续模型的方法。包括：1、插值方法 2、数据拟合方法 3、最佳逼近要求：1、了解MATLAB功能函数 2、掌握数据模拟的各种方法及应用,3.1问题的提出,例一：在某化学反应里，侧的生成物的质量浓度y与时间t(min)的关系入表。为了研究该化学反应的性质，如反映速率等，欲求y与t之间的关系y=f(t)。,例二：逢山开道（山区修建公路的路线选择问题1994）在某山区修建公路，已知该山区的地形高度见表一（单位m）。雨季在山谷形成一溪流，雨量最大时溪流水面宽度W与x（溪流最深处的）坐标的关系可以表示为,要求：从山脚开始经居民点至矿区修一条公路（一般公路、桥、隧道）给出路线设计方案使工程成本最小。,部分数据表：,2000,工程要求：,例三：水道测量模型（1989.美国）问题 水平方向的坐标x，y以Td（=0.914m）为单位,水深方向Z 以Ft(=30.48cm)为单位.下表给出了水面直角坐标(x,y)处的水深Z，这是在低潮时测得的。如果船的吃水深度为5Ft，试问在矩形城(75，200)(-50，50)中行船应避免进入那些区域?,工作：1、要根据数据做出地貌图（三维）,等高线（等势图）二维。2、由坡度限制，沿给定的网格点选择路线,工作：将海底曲面图绘制出来。,总之：根据给定采样数据而定出其实际分布图，这就是数值模拟问题，采用的方法就是插值和拟合。,3.2 数值模拟概论,在科学与工程等实际问题中，实际数学模型或者难于建立或者难于求解，但其数据模型（由实验或测量所得到的一批离散数据）容易得到。能否通过处理这些数据来建立实际模型呢？这里我们仅以一维问题来说明。,给定：y=f(x)数据,通过这批数据希望找到对象y=f(x)的更多的信息（或全部信息）。通常只能找到f(x)的近似表达式(x),(x)是根据研究对象的特征确定的。,对象是指数增（减）并最终稳定到某个值,则,对象是直线运动、匀加速运动，则,一般地，我们都要通过对研究对象实际的了解，做以下工作：,1、确定它的基本特征函数,2、对这些特征函数（基函数）进行一种恰当的构造，得到f(x)近似的数学模型:,最常用也就是最简单的一种构造方式就是：,3、按某种寻优策略，由已知数据确定未知参数,寻优策略的不同，产生了不同的数值逼近方法,3.3 插值方法,寻优策略就是过点，即,根据需要可附加补充原则:,1、p阶光滑度,2、边界条件：,处光滑性,周期性等,（整体）误差：,插值方法中最常用的一类就是代数插值。,代数插值：即多项式的插值,一、代数插值,（一）拉格朗日插值公式（Lagrange),给定：y=f(x)数据,确定次数不超过n的多项式P(x),使,拉格朗日插值公式的讨论,线性插值公式（n=1的情况）,显然，两插值节点越近，越精确。,2、二次插值（抛物线插值，n=2的情况）,n越大精度越高！是这样吗？下面我们看一实例。,实例演示：分别用n=1,2,4,8,16,32时的拉格朗日插值逼近f(x)。,function f=lagelangri(t,y,x)n=length(t);for j=1:n l(j)=1;for i=1:n if i=j l(j)=l(j)*(x-t(i)/(t(j)-t(i);end endendf=y*l;,定义lagrange函数,function f=shiyanhs(x)f=1./(1+25*x.2);,定义实验函数,function g=lagelangritu(a,b,n)h=(b-a)/n;t=a:h:b;f0=shiyanhs(t);x=a:0.02:b;m=length(x);for i=1:m p(i)=lagelangri(t,f0,x(i);end f=shiyanhs(x);plot(x,f,r,x,p);gtext(f(x);gtext(p(x,n);axis(-1,1,-1,1);,画图：f(x)和p(x),总之，我们看到利用拉格朗日公式逼近函数f(x),n小不行，n大也不行。这是为什么呢？,综合以上分析我们可以得到如下结论：在小区间上应用低次拉格朗日差值公式（n=6)有效。,问题是对大的区间怎么办？那就是用分段低次插值。,（二）分段低次插值,分段线性插值:,给定：y=f(x)数据,确定P(x),使,根据要求立即得到：,若用基函数的形式表示，则：,function f=fenduanxianxing(t,y,x)%x在t(1),t(n)n=length(t);i=1;while 1 if(x=t(i),定义分段线性插值函数,实验：用分段线性插值逼近函数,function g=fenduanxianxingtu(a,b,n)h=(b-a)/n;t=a:h:b;f0=shiyanhs(t);x=a:0.02:b;m=length(x);for i=1:m f(i)=fenduanxianxing(t,f0,x(i);end f1=shiyanhs(x);plot(x,f1,.-r,x,f);gtext(f(x);gtext(p(x,n);axis(-1,1,-1,1);,g=fenduanxianxingtu(-1,1,50),计算：,分段线性插值的讨论：,误差分析,下面的实例演示可以说明,g=fenduanxianxingtu(-1,1,n),分别用n=4，8,16,32的分段线性差值逼近函数（h=2/n),样条函数与样条插值：,给定：y=f(x)数据,确定s(x),使,若s(x)存在，则称s(x)为对应于给定数据的一个k次样条插值。满足条件1和2的函数S(x)称为k次样条函数。,样条函数的这种数学定义来源于 绘图员沿着受压铁约束的样条（一根具有很好弹性的木条），画出各种光滑曲线（样条曲线）。这种样条曲线可以看成弹性细梁受集中载荷作用而生成的挠度曲线。,S(x)的存在性讨论：,实际中常用k=2,3的情况，即二次样条和三次样条,二次样条插值公式,因为二次样条函数的确定需要n+2个条件，而,给定了n+1个，所以还需补充一个条件（自然是边界条件）。通常有,补充条件1下的二次样条插值公式,采用如下方法来确定：,联立式（1）、（2）、（3）确定了二次样条插值公式：,这样不如前面的方法简单！,实验：用二次样条插值逼近f(x),h=(b-a)/n.比较n=4,8,16,32的效果,二次样条插值的误差分析,三次样条插值公式,因为三次样条函数的确定需要n+3个条件，而,给定了n+1个，所以还需补充二个条件（自然是边界条件）。通常有,解由插值条件和补充条件构成的线性方程组，可以得到三次样条插值公式。但这样比较复杂！下面只就问题2来讨论，并采用另一种途径称之为三弯矩法。,综合上述，构造三次样条插值公式的算法为：第一步：式（2）、（3）构成如下方程组,三次样条插值的误差分析,实验：用三次样条插值逼近f(x),h=(b-a)/n.比较n=4,8,16,32的效果,均匀分划下的样条函数与样条插值,一、基本B样条函数从单位阶跃函数谈起,0阶B样条函数：,1阶B样条函数：,2阶B样条函数：,3阶B样条函数：,一般的有k次B样条函数：,实际上常用k=1,2,3。而k=3最用实用价值。,二、均匀分划下的样条函数逼近,已知,B样条函数插值问题：,问题的解：,k=1时，一次B样条插值公式：,它就是分段线性插值。,k=2时，二次B样条插值公式：,解,得到,误差分析结果与前面一般二次样条插值一样！但在边界点的光滑性不能保证。,k=3时，三次B样条插值公式：,解,得到,误差分析结果与前面一般三次样条插值一样！但在边界点的光滑性不能保证。,注：,为使B样条插值在边界点也光滑，可以在两边延拓。以三次B样条插值为例。即,解：,实验：用一、二、三阶B样条插值逼近f(x),h=(b-a)/n.比较n=4,8,16,32的效果,n为4的情况,n=8的情况,n=16的情况,n=32的情况,阶梯函数与折线函数的磨光,定理1：,磨光法在外形设计问题重的应用,例 某外形设计问题所给的型值 如下表。要求寻找这样的曲线，满足：1）通过两个端点，2）在x(i)处逼近原型值的绝对误差不超过1厘米，3）在x(i)处的二阶导数符号与原型值点的二阶差符号相同。,下面我们采用盈亏修正法来提高精度。其算法如下：,实例计算结果,3.4 数据拟合的最小二乘法,问题的引入,先讨论两个变量x,y的情况。通过观测变量x,y积累了一组资料（x(i),y(i),i=1:n，一般来说n比较大。我们的任务是从这批实验数据出发，寻求一近似函数(x)曲逼近y。由于观测数据都带有观测误差，数目又较大，对于这类问题运用插值函数取描述y往往是不适当的。现在，用开始引入的例1来说明建立近似函数的一种方法，即最小二乘法。,研究y与t的关系的最小二乘法，通常步骤如下：,1）先描t,y的草图,2）凭视觉猜想一数学模型,3）计算总误差,4）希望总误差最小,数据拟合的最小二乘法,对称矩阵方程,例1的数据拟合结果,显然，比第一种模型要好得多！,