基于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,、当,0eyta,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;,end,lmta=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=lmta,end,动态位移反幂法,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);,end,end,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-1k: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);,end,for,i=1:n-1,for,j=i+1:n,r(i,j)=a(i,j);,end,end,for,i=1:n-1,r(i,i)=-arfa(i);,end,r(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;,end,h=a;,for,i=1:n-2,h(i+1,i)=-arfa(i);,end,for,j=1:n-2,for,i=j+2:n,h(i,j)=0;,end,end,注：当,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,ss,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),。,1 2 3 4 6 8 10 12 14 16,4.00 6.41 8.01 8.79 9.53 9.86 10.33 10.42 0.53 10.61,例二：,逢山开道（山区修建公路的路线选择问题,1994,）,在某山区修建公路，已知该山区的地形高度见表一（单位,m,）。雨季在山谷形成一溪流，雨量最大时溪流水面宽度,W,与,x,（溪流最深处的）坐标的关系可以表示为,要求：从山脚开始经居民点至矿区修一条公路（一般公路、桥、隧道）给出路线设计方案使工程成本最小。,3200,1430,1450,1460,1550,1500,1600,2800,950,1190,1370,1500,1200,1100,2400,910,1090,1270,1500,1200,1100,2000,880,1060,1230,1390,1500,1500,1600,830,980,1180,1320,1450,1420,1200,740,880,1080,1130,1250,1280,800,*650,760,880,970,1020,1050,400,510,620,730,800,850,870,0,370,470,550,660,670,690,Y/x,0,400,800,1200,1600,2000,部分数据表：,工程种类,一般路段,桥梁,隧道,工程成本（元,/m,）,300,1500,（长度,300m,）,=,上升高度,/,水平距离,=0.125,=0,=0.100,2000,工程要求：,例三,：水道测量模型（,1989.,美国）,问题 水平方向的坐标,x,，,y,以,Td,（,=0.914m,）为单位,水深方向,Z,以,Ft(=30.48cm),为单位,.,下表给出了水面直角坐标,(x,y),处的水深,Z,，这是在低潮时测得的。如果船的吃水深度为,5Ft,，试问在矩形城,(75,，,200)(-50,，,50),中行船应避免进入那些区域,?,工作：,1,、要根据数据做出地貌图（三维）,等高线（等势图）二维。,2,、由坡度限制，沿给定的网格点选择路线,x,129.0,140.0,108.5,88.0,185.5,195.5,105.5,y,7.5,141.5,28.0,147.0,22.5,137.5,85.5,z,4,8,6,8,6,8,8,157.5,107.5,77.0,81.0,162.0,162.0,117.5,-6.5,-81.0,3.0,56.5,-66.5,84.0,-38.5,9,9,8,8,9,4,9,x,y,z,工作：将海底曲面图绘制出来。,总之：根据给定采样数据而定出其实际分布图，这就是数值模拟问题，采用的方法就是插值和拟合。,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,end,end,f=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=t(i)&(x=t(i+1),break,end,;,i=i+1;,end,f=yi)*(x-t(i+1)/(t(i)-t(i+1)+y(i+1)*(x-t(i)/(t(i+1)-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(-11,-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),处的二阶导数符号与原型值点的二阶差符号相同。,X(i)cm,30 80 130 180 230 280 330 380 430 480,Y(i)cm,80 110 132 148.75 163 175 185.5 195 204 212.75,下面我们采用盈亏修正法来提高精度。其算法如下：,实例计算结果,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,的数据拟合结果,显然，比第一种模型要好得多！,