数值分析实验教程

TOC\o"1-3"\h\z第二篇数值分析 3第1章 绪论 41.1 2的平方根计算 41.2 计算效率探讨 5实验题 7第2章 插值法 82.1 拉格朗日插值多项式的存在性 82.2 利用拉格朗日插值多项式计算函数值 82.3 差商表构造 102.4 利用牛顿插值多项式计算函数值 112.5 龙格现象 122.6 分段线性插值的逼近性 142.7 拉格朗日插值多项式与埃米特插值多项式的比较 162.8 拉格朗日插值多项式与三次样条插值函数的比较 18实验题 22第3章 函数的数值逼近 243.1 伯恩斯坦多项式逼近连续函数的动画演示 243.2 函数的最佳平方逼近多项式 263.3 希尔伯特矩阵的病态性 273.4 多项式拟合模型的选取 283.5 非线性拟合模型的选取 29实验题 31第4章 数值积分 324.1 牛顿—科茨公式的数值稳定性 324.2 复化求积思想的动画演示 344.3 自动选取步长的复化梯形求积算法 364.4 龙贝格数表 374.5 龙贝格求积方法 39实验题 41第5章 非线性方程求根 435.1 两分法求方程的根 435.2 定积分中值定理的几何证明 455.3 迭代法的敛散性研究 475.4 艾特肯方法 495.5 牛顿法 505.6 弦截法 515.7 快速弦截法 53实验题 54第6章 常微分方程数值解 566.1 欧拉方法的几何意义 566.2 欧拉方法与改进的欧拉方法的比较 576.3 四阶经典公式 596.4 四阶阿当姆斯预报校正系统 616.5 变步长的龙格库塔方法 63实验题 66第7章 线性方程组求解 687.1 高斯顺序消去法 687.2 高斯列主元消去法 707.3 高斯全主元消去法 727.4 方阵的LU分解 767.5 选方阵列主元的LU分解 787.6 追赶法 807.7 平方根法 827.8 改进的平方根法 847.9 雅可比迭代 857.10 高斯—塞德尔迭代 88实验题 89

第二篇数值分析在信息时代，科学技术出现了前所未有的发展，其中数学应用的广泛性和深入性已成为现代科技发展的重要特征。数学与科学计算、理论研究、科学实验并列为科学研究的三大支柱，计算机已成为科学与工程技术等领域不可缺少的工具。因此，在大学数学教育中，应特别注重对实际问题的数学建模，为数学模型求解提供有效的计算算法以及运用计算机进行求解等能力的培养与训练。数值分析是研究各种数学问题求数值解的方法，离散化、递推化是它处理问题的主要手段，误差分析是它研究的核心问题，以计算机和数学软件为工具进行数值计算是它的显著特征。MATLAB具有强大的数值计算功能，可以进行数值插值、曲线拟合、数值积分、非线性方程求根、线性方程组求解、微分方程求数值解等运算。本篇主要给出了数值分析涉及到数学实验，这些实验的编程思想和MATLAB程序代码。通过本篇的学习，能对数值分析中主要概念、主要数值算法有更进一步的认识，更好地掌握MATLAB编程，并学会利用MATLAB解决科学计算的问题。

绪论2的平方根计算实验用例用递推算法求的近似值，并回答下述问题1.的有效数字是几位；2.的有效数字的位数是否为的2倍？若是，这是否表明该算法具有较高的计算效率；3.若近似值的误差精度要求为，编程计算的近似值。实验步骤及说明1.在MATLAB命令空间键入formatlong，按下Enter键，设置数值格式为实长型。2.键入sqrt(2)并按下Enter，可得到具有较高精度的近似值，我们假定该值为的精确值，其值为1.41421356237310。3.键入x=2，按下Enter键，再键入x=(x+2/x)/2，再按下Enter键，得到第1次递推值x1；按下↑键，调出表达式x=(x+2/x)/2，按下Enter键，得到第2次递推值x2；再按下↑键，调出表达式x=(x+2/x)/2，再按下Enter键，得到第3次递推值x3。将该值与的精确值作比较，可知x具有6位有效数字。4.继续按下↑键，调出表达式x=(x+2/x)/2，并按下Enter键，得到第4次递推值x4。将该值与的精确值进行比较，可知x4具有12位有效数字。这表明，第4次递推值的有效数字的位数比第3次递推值的有效数字的位数增加了2倍，该算法确实具有较高的计算效率。为了解决实验中的第3个问题，我们需要编写MATLAB程序。编程思想1.从数学上可以证明，序列单调下降趋向于。与的误差，可以用来近似代替。一旦，便终止递推过程，取作为的近似值。为了方便，我们引入变量err来保存误差。2.由于递推过程仅需要保留与两个值，我们分别用x1与x2来代替。3.递推开始时，我们设置误差err的初值为1。事实上，只要是满足条件err>ε的任意err值都可以作为初值。程序代码eps=10^(-8);%误差精度err=1;%初始误差x1=2;%递推初值whileerr>=epsx2=(x1+2/x1)/2;%递推err=abs(x2-x1);%误差计算x1=x2;%准备下一次的递推初值endfprintf('2的平方根的近似值为%10.8f\n',x1)%输出近似值运行结果2的平方根的近似值为1.41421356计算效率探讨实验用例利用函数展开式可得到计算ln2的级数利用函数展开式可得到计算ln2的另一个级数1.给出（1）式前15项的部分和，（2）式前3项的部分和，指出计算结果作为的近似值，各含几位有效数字，哪一个更好；2.再计算（2）式前4项，前5项的值，每增加求和表达式中的一项，其近似值的有效数字就增加一位；3.用级数计算近似值收敛缓慢的原因是，求和算式中出现了正负项交错，使得部分和数列在附近左右跳动，试通过编程验证这一现象。实验步骤及说明1.在MATLAB命令空间键入formatlong，按下Enter键，设置数值格式为实长型。键入log(2)并按下Enter，可得到ln2的近似值约为0.69314718055995。2.键入s1=1-1/2+1/3-1/4+1/5-1/6+1/7-1/8+1/9-1/10+1/11-1/12+1/13-1/14+1/15，按下Enter键，得到s1=0.72537185037185，该值具有0位有效数字。3.键入s2=(2/3)*(1+(1/3)^2/3+(1/3)^4/5)，按下Enter键，得到s2=0.69300411522634，该值具有3位有效数字。由此可见，（2）式要比（1）更好。4.键入s3=(2/3)*(1+(1/3)^2/3+(1/3)^4/5+(1/3)^6/7)，得到s3=0.69313475733229，该值具有4位有效数字。s4=(2/3)*(1+(1/3)^2/3+(1/3)^4/5+(1/3)^6/7+(1/3)^8/9)，得到s4=0.69314604739083,该值具有5位有效数字。编程思想1.输入部分和的项数n，并将部分和存入数组s1,s2,…,sn，其部分和的算法如下：（1）t←1，s1←1（2）对于k=2,3,…,n，做以下操作①t←(-1)t②sk←sk-1+t/k（3）输出sn2.作数组s1,s2,…,sn以及y=ln2的图像。程序代码close%关闭所有图形窗口clearall%清除命令空间中所有变量n=input('请输入部分和的项数n=');t=1;%符号变量初始化s(1)=1;%部分和数列第1项fork=2:1:nt=(-1)*t;%符号交错s(k)=s(k-1)+t/k;%部分和计算endfprintf('计算得到的近似值s(%4d)=%10.9f\n',n,s(n))%输出近似值holdonk=1:1:n;plot(k,s,'b')%作部分和数列图像plot([1,n],[log(2),log(2)],'r')%作y=ln2的图像holdoff运行结果请输入部分和数列的项数n=100计算得到的近似值s(100)=0.688172179部分和数列的折线图如图1-1所示。图1-1部分和数列折线图实验题1.编写秦九韶算法程序，并用该程序计算多项式在自变量时的值。2.考虑无穷级数，它是微积分中著名的发散级数。在计算上计算该级数的部分和，会得到怎样的结果？为什么？3.对于，用级数展开式近似计算函数值，这样的算法结果是否会好？为什么？4.考虑数列，它的统计平均值为,它的标准差为，它等价于，作为标准差的两种算法，从数值计算的角度，你对它们作何种评价。

插值法拉格朗日插值多项式的存在性实验用例利用求解线性方程组的方法，给出适合插值条件xi－1125yi－77－435的三次插值多项式。实验步骤及说明1.在MATLAB命令空间依次键入以下命令formatrat%设置数据为有理数类型x=[-1,1,2,5]';y=[-7,7,-4,35]';A=[x.^3，x.^2，x.^1,x.^0];A\y可以得到ans=[2，-10，5，10]'。这表明，满足上述插值条件的拉格朗日插值多项式为2.在输入了插值条件列向量x，y之后，调用函数v=vander(x)，可得到相应的范德蒙矩阵v=-11-11111184211252551再输入命令v\y，也可得到这一结果。3.在输入了插值条件向量x，y之后，运用MALTAB内建函数polyfit，键入命令polyfit(x,y,3)也可很方便地获得这一结果。利用拉格朗日插值多项式计算函数值实验用例求满足插值条件xi100121144yi101112的拉格朗日插值多项式在x=115处的值。实验步骤及说明1.利用MATLAB函数polyfit，polyval可给出拉格朗日插值多项式，并计算该多项式在某点的函数。解决上述实验问题，只需在命令空间依次键入以下命令x=[100,121,144];y=[10,11,12];polyval(polyfit(x,y,2),115)可以得到ans=10.7228。2.利用拉格朗日插值基函数表示形式编程计算，其算法如下：(1)将插值节点存入向量x，其对应的函数值存入向量y(2)测出x所含节点个数n(3)输入z(4)L←0(5)对于i=1,2,…,n，做以下操作①t←1②对于j=1,2,…,n，做以下操作(ⅰ)若j≠i，则t←t(z-xj)/(xi-xj)③L←L+tyi(6)输出L程序代码x=[100,121,144];y=[10,11,12];n=length(x);z=input('请输入欲求其函数值的节点z=');L=0;fori=1:nt=1;forj=1:nifj~=it=t*(z-x(j))/(x(i)-x(j));%计算拉格朗日插值基函数endendL=L+t*y(i);%计算拉格朗日插值多项式endfprintf('函数值L(%8.4f)=%10.6f\n',z,L)运行结果请输入欲求其函数值的节点z=115函数值L(115.0000)=10.722756差商表构造实验用例对于插值条件xi0.40.550.650.800.901.05yi0.410750.578150.696750.888111.026521.25382编程构造差商表。编程思想1.将插值节点存入向量x，其对应的函数值存入向量y。2.n←length(x)3.newton←[x',y']4.对于j=2,…,n，做以下操作(1)对于i=n,n-1,…,2,1，做以下操作①若i≥j，则yi←(yi–yi-1)/(xi–xi-j+1)；否则，yi←0(2)newton←[newton,y']5.显示牛顿差商表newton。需要说明的是，差商表上三角部分的0仅表示空格，不具有任何实际意义。向量y中最后存放的实际上是0阶，1阶，…，n-1阶差商。程序代码x=[0.4,0.55,0.65,0.80,0.90,1.05];%插值节点y=[0.41075,0.57815,0.69675,0.88811,1.02652,1.25382];%插值节点处的函数值n=length(x);%测x所含节点个数newton=[x',y'];%给出牛顿差商表的前两列forj=2:n%计算第j-1阶差商fori=n:-1:1ifi>=jy(i)=(y(i)-y(i-1))/(x(i)-x(i-j+1));elsey(i)=0;%0不具有实际意义,仅表示空格endendnewton=[newton,y'];%将j-1阶差商存入牛顿差商表的第j+1列enddisp('下三角状的牛顿差商表如下:')disp(newton)运行结果z=0.40000.4108000000.55000.57821.116000000.65000.69671.18600.28000000.80000.88811.27570.35890.1973000.90001.02651.38410.43350.21300.031201.05001.25381.51530.52490.22870.03140.0003利用牛顿插值多项式计算函数值实验用例已知函数的部分函数值xi0.40.550.650.800.901.05yi0.410750.578150.696750.888111.026521.25382利用牛顿插值多项式，计算的近似值并估计其误差。编程思想1.将插值节点存入向量x，其对应的函数值存入向量y，z←0.596。2.将1阶，2阶，…，n-1阶差商依次存入向量y中的y2，y3，…，yn，其算法如下：(1)对于j=2,3,…,n，做①①对于i=n,n-1,…,j，做yi←(yi–yi-1)/(xi–xi-j+1)3.用u表示牛顿插值多项式在z处的值，u的计算采用秦九韶算法。(1)u←yn(2)对于i=n-1,…,2,1，做u←yi+u(z-xi)4.用r表示误差，r计算的算法如下。(1)r←yn(2)对于i=1,2,…,n,做r←r(z-xi)5.输出u，r。程序代码x=[0.4,0.55,0.65,0.80,0.90,1.05];y=[0.41075,0.57815,0.69675,0.88811,1.02652,1.25382];z=0.596;%给出需求其函数值的点n=length(x);forj=2:1:nfori=n:-1:jy(i)=(y(i)-y(i-1))/(x(i)-x(i-j+1));%计算差商endendu=y(n);fori=n-1:-1:1u=y(i)+u*(z-x(i));%计算牛顿插值多项式的值endr=y(n);fori=1:nr=r*(z-x(i));%误差估计endfprintf('f(%8.5f)=%8.5f\n',z,u);%输出函数近似值fprintf('截断误差约为R(%8.5f)=%15.14f\n',z,abs(r));%输出近似误差运行结果f(0.59600)=0.63192截断误差R(0.59600)=0.00000000401693龙格现象实验用例对于被插函数以及插值节点1.作出与在[-5,5]上的图像，观察龙格现象；2.如果被插函数为或，试利用图像来说明龙格现象是否会发生。编程思想1.为被插函数创建函数文件f.m，其代码如下：functiony=f(x)y=1./(1+x.^2);2.为使程序具有一般性，用实型变量a，b分别表示插值区间的左右端点，这里，a=-5，b=5。用正整数变量n表示区间[a,b]的等分数，这里n=10，而插值节点xi=a+(b-a)i/n(i=0,1,…,n)。3.计算y=Ln(x)在xi=a+(b-a)i/(4n)(i=0,1,…,4n)处的函数值，其程序的核心代码参见实验2.2。4.分别用蓝色和红色线分别作函数y=1/(1+x2)和y=Ln(x)的图像。程序代码closeclearalla=input('请输入插值区间的左端点a=');b=input('请输入插值区间的右端点b=');n=input('请输入插值区间的等分数n=');x=a:(b-a)/n:b;%给出插值条件y=f(x);z=a:(b-a)/(4*n):b;%给出需计算其插值多项式值的节点m=length(z);fork=1:m%计算插值多项式值L(k)=0;fori=1:n+1t=1;forj=1:n+1ifj~=it=t*(z(k)-x(j))/(x(i)-x(j));endendL(k)=L(k)+t*y(i);endendplot(x,y,'b')%作被插函数图像holdonplot(z,L,'r')%作拉格朗日插值多项式图像title('龙格现象演示')%给图形窗口加标题holdoff运行结果请输入插值区间的左端点a=-5请输入插值区间的右端点b=5请输入插值区间的等分数n=10程序运行后的图像见图2-1。图2-1龙格现象的演示实验第2部分只需将函数文件f.m的内容分别修改为functiony=f(x)y=exp(-abs(x));或functiony=f(x)y=sin(x);输入内容仍为a=-5，b=5，n=10可得到以下图像。观察图像可发现，对于被插函数，仍有龙格现象出现。而对于被插函数，不发生龙格现象。程序运行结果参见图2-2。图2-2检验龙格现象是否发生分段线性插值的逼近性实验用例对于被插函数以及插值节点1.当用户输入插值区间的等分数n之后，作出被插函数与它的分段线性插值函数的图像；2.指出当n增大时，是否会发生龙格现象；3.观察n=2，4，8，16，32时的图像，你可以得出怎样的结论。编程思想1.为被插函数创建函数文件f.m，其代码如下：functiony=f(x)y=1./(1+x.^2);2.程序的主要结构与实验2.5基本相同，主要区别仅在分段线性插值函数y=S1(x)在节点xi=a+(b-a)i/(2n)(i=0,1,…,2n)处的函数值计算，其主要算法如下：(1)取z=a:(b-a)/(2n):b(2)测向量z的长度并赋给m(3)对于j=1,2,…,m，做以下操作①对于i=1,2,…,n，做以下操作(ⅰ)若xi≤zj≤xi+1S1j←yi(zj–xi+1)/(xi–xi+1)+yi+1(zj–xi)/(xi+1-xi)break程序代码closeclearalla=input('请输入插值区间的左端点a=');b=input('请输入插值区间的右端点b=');n=input('请输入插值区间的等分数n=');x=a:(b-a)/n:b;%插值节点y=f(x);%插值节点处的函数值z=a:(b-a)/(2*n):b;%给出需计算其分段线性插值函数值的节点m=length(z);forj=1:m%分段线性插值函数在节点z处的函数值计算fori=1:nifz(j)>=x(i)&z(j)<=x(i+1)S1(j)=y(i)*(z(j)-x(i+1))/(x(i)-x(i+1))+y(i+1)*(z(j)-x(i))/(x(i+1)-x(i));breakendendendplot(z,S1,'r')%作分段线性插值函数图像holdonu=a:0.01:b;v=f(u);plot(u,v,'b')%作被插函数图像title('分段线性插值函数逼近性演示')holdoff运行结果请输入插值区间的左端点a=-5请输入插值区间的右端点b=5请输入插值区间的等分数n=10程序运行后的图像见图2-3。图2-3分段线性插值函数逼近性演示反复运行该程序，可以发现：当n增大时，龙格现象不会发生；当n=2,4,8,16,32时，分段线性插值函数的图像逐渐逼近被插函数的图像。拉格朗日插值多项式与埃米特插值多项式的比较实验用例给定函数，节点以及插值条件1.作出被插函数，10次拉格朗日插值多项式，21次的埃米特插值多项式的图像；2.观察图像，你认为结论“在插值区间的某些点处，21次埃米特插值多项式比10次拉格朗日插值多项式的龙格现象更加显著”是否正确。编程思想1.将插值节点存入向量x，将其对应的函数值、导数值分别存入向量y、ybar，n=length(x)。2.给出需计算拉格朗日多项式，埃米特插值多项式值的节点z=-5:0.01:5，m=length(z)。3.对于k=1,2,…,m，计算拉格朗日插值多项式的值L(zk)，埃米特插值多项式的值H(zk)。4.作被插函数，拉格朗日插值多项式，埃米特插值多项式的图像。下面对L(zk)、H(zk)计算的原理与编程思想作说明。数值分析教材中给出了拉格朗日插值多项式与埃米特插值多项式的计算公式其中，。其中，，对这两个计算公式，可按以下算法来编程计算。(1)对于zk（k=1,2,…,m），做以下操作①Lk←0，Hk←0②对于i=1,2,…,n，做以下操作(ⅰ)t←1，s←1(ⅱ)对于j=1,2,…,n，做以下操作若j≠i，则t←t(zk-xj)/(xi-xj)，s←s+2(xi-zk)/(xi-xj)(ⅲ)Lk←Lk+yit，Hk←Hk+yist2+yi’(zk–xi)t2程序代码closeclearall%给出插值条件x=-5:1:5;y=1./(1+x.^2);ybar=-2*x./(1+x.^2).^2;n=length(x);%计算拉格朗日插值多项式，埃米特插值多项式在点z=-5:0.01:5的值z=-5:0.01:5;m=length(z);fork=1:mL(k)=0;%拉格朗日插值多项式赋初值H(k)=0;%埃米特插值多项式值赋初值fori=1:nt=1;%拉格朗日插值基函数赋初值s=1;%埃米特插值基函数Alfa中的和式赋初值forj=1:nifj~=it=t*(z(k)-x(j))/(x(i)-x(j));s=s+2*(x(i)-z(k))/(x(i)-x(j));endendL(k)=L(k)+y(i)*t;H(k)=H(k)+y(i)*s*t^2+ybar(i)*(z(k)-x(i))*t^2;endendu=-5:0.01:5;v=1./(1+u.^2);plot(u,v,'b')%作被插函数图像holdonplot(z,L,'r')%作拉格朗日插值多项式图像plot(z,H,'k')%作埃米特插值多项式图像holdoff运行结果图2-4拉格朗日插值多项式与埃米特插值多项式的比较图中蓝色线表示被值函数，红色线表示拉格朗日插值多项式，黑色线表示埃米特插值多项式。从图中可发现，埃米特插值多项式在0的附近与被插函数吻合的较好，而在-5，5附近会出现更为显著的龙格现象。拉格朗日插值多项式与三次样条插值函数的比较实验用例给定函数，节点以及插值条件1.计算处的拉格朗日插值多项式，三次样条插值函数的值，并对计算进行比较；2.作被插函数，拉格朗日插值多项式，三次样条插值函数的图像。观察图像，你认为“三次样条插值函数比拉格朗日插值多项式能更好地逼近被插函数”这一结论是否正确。编程思想1.将插值节点存入向量x，其对应的函数值存入向量y，n=length(x)。将两个端点的导数值存入向量ybar。取z=-5:0.5:5，m=length(z)。2.计算L(x)在z处的函数值。3.计算S(x)在z处的函数值。4.输出节点，被插函数值，三次样条插值函数值以及拉格朗日插值多项式的值。5.作被插函数，拉格朗日插值多项式，三次样条插值函数的图像。L(x)在z处的函数值计算，可参考实验2.2或实验2.7来编程。下面着重介绍S(x)在z处函数值计算的原理以及算法。数值分析教材中给出的第一边界条件下三次样条插值函数计算公式如下：其中，矩满足以下三弯矩方程这里，，，，，，其计算的算法如下：(1)λ1←1，μn←1(2)对于i=1,2,…,n–1，做以下操作①hi←xi+1–xi(3)d1←6((y2–y1)/h1–ybar1)/h1，dn←6(ybar2–(yn–yn-1)/hn-1)/hn-1(4)对于i=2,3,…,n-1，做以下操作①λi←hi/(hi-1+hi)②μi←hi-1/(hi-1+hi)③di←6((yi+1–yi)/hi–(yi–yi-1)/hi-1)/(hi-1+hi)(4)用追赶法求解三弯矩方程，计算矩M1,M2,…,Mn（这一算法参见实验7.4）(5)对于zk，k=1,2,…,m，做以下操作①对于i=1,2,…,n–1，做以下操作(ⅰ)若xi≤zk≤xi+1，则Sk←Mi(xi+1–zk)3/(6hi)+Mi+1(zk–xi)3/(6hi)+(yi–Mihi2/6)(xi+1–zk)/hi+(yi+1–Mi+1hi2/6)(zk–xi)/hibreak程序代码closeclearall%给出插值条件x=-5:1:5;y=1./(1+x.^2);n=length(x);ybar=[-10/(26)^2,10/(26)^2];%给出需计算函数值的节点z=-5:0.5:5;m=length(z);%计算L(x)在z处的函数值fork=1:mL(k)=0;fori=1:nt=1;forj=1:nifj~=it=t*(z(k)-x(j))/(x(i)-x(j));endendL(k)=L(k)+y(i)*t;endend%计算S(x)在z处的函数值%计算三弯矩方程的系数lmd(1)=1;mu(n)=1;fori=1:n-1h(i)=x(i+1)-x(i);endd(1)=6*((y(2)-y(1))/h(1)-ybar(1))/h(1);d(n)=6*(ybar(2)-(y(n)-y(n-1))/h(n-1))/h(n-1);fori=2:n-1lmd(i)=h(i)/(h(i-1)+h(i));mu(i)=h(i-1)/(h(i-1)+h(i));d(i)=6*((y(i+1)-y(i))/h(i)-(y(i)-y(i-1))/h(i-1))/(h(i-1)+h(i));end%追赶法解三弯矩方程，求矩Mbata(1)=lmd(1)/2;fori=2:n-1bata(i)=lmd(i)/(2-mu(i)*bata(i-1));endu(1)=d(1)/2;fori=2:nu(i)=(d(i)-mu(i)*u(i-1))/(2-mu(i)*bata(i-1));endM(n)=u(n);fori=n-1:-1:1M(i)=u(i)-bata(i)*M(i+1);end%计算三次样条插值函数的值fork=1:mfori=1:n-1ifz(k)>=x(i)&z(k)<=x(i+1)S(k)=M(i)*(x(i+1)-z(k))^3/(6*h(i))+M(i+1)*(z(k)-x(i))^3/(6*h(i))+(y(i)-M(i)*h(i)^2/6)*(x(i+1)-z(k))/h(i)+(y(i+1)-M(i+1)*h(i)^2/6)*(z(k)-x(i))/h(i);breakendendend%计算f(x)在z处的函数值F=1./(1+z.^2);%输出节点，函数值，样条值以及拉氏插值disp('节点函数值样条值拉氏插值')disp([z',F',S',L'])%作图plot(z,F,'b')%用蓝线作被插函数图像holdonplot(z,L,'k')%用黑线作拉格朗日插值多项式图像plot(z,S,'r')%用红线作三次样条插值函数图像holdoff运行结果运作程序，可得到以下输出值以及图2-5。由于被插函数是偶函数，我们仅给出了z取负值时所对应的函数值。节点函数值样条值拉氏插值-5.00000.03850.03850.0385-4.50000.04710.04721.5787-4.00000.05880.05880.0588-3.50000.07550.0748-0.2262-3.00000.10000.10000.1000-2.50000.13790.14000.2538-2.00000.20000.20000.2000-1.50000.30770.29740.2353-1.00000.50000.50000.5000-0.50000.80000.82050.843401.00001.00001.0000图2-5拉格朗日插值多项式与三次样条插值函数的比较从数值计算结果与图像可以看出，三次样条插值函数S(x)比拉格朗日插值多项式L(x)能更好地逼近被插函数f(x)。实验题1.设，，而可取0.5，0.6，0.7中的任一点，根据下表x0.40.50.60.70.8sinx0.38940.47940.56460.64420.7174使用二次拉格朗日插值函数求的近似值时，问(1)如何选节点可使其近似值的误差较小(要求给出理论证明)；(2)通过数值实验，完成下述表格，并验证你的理论结果。x10.50.60.7L2(0.6389)|R(0.6389)|2.已知的函数表x0.500.550.600.650.70y0.106530.02695-0.05119-0.12795-0.20342(1)作函数在(0，1)上的图像，说明方程是否有根。若有根，给出它的大致区间，并运用严格的数学方法证明该方程的根存在且唯一；(2)运用牛顿插值，给出方程位于区间(0，1)内根的近似值，并估计误差。3.已知函数表如下，利用三次样条插值函数，求该函数在的近似值。x0123y0236y’104.根据以下函数表，求三次样条插值函数。x27.7282930y4.14.34.13.0y’3.0-4.05.编程实现分段线性插值函数逼近被插函数的动画。6.区间[a,b]上的切比雪夫点定义为对于被插函数比较以切比雪夫点、等距点为插值节点的拉格朗日插值多项式的图像。

函数的数值逼近伯恩斯坦多项式逼近连续函数的动画演示实验用例伯恩斯坦多项式其中，用动画形象地证明，在[0，1]上，一致逼近函数（可取或其它函数来进行实验）。编程思想1.为被逼近函数创建函数文件f.m，其代码如下：functiony=f(x)y=exp(x);2.程序的主要算法如下：（1）x=0:0.05:1，y=f(x)，q←length(x)（2）输入伯恩斯坦多项式的次数数组n，m←length(n)（3）对于i=1,2,…,m，做以下操作①作被逼近函数的图像，并固化图形②对于j=1,2,…,q，计算伯恩斯坦多项式在xj处的值Bj③作伯恩斯坦多项式的图像④图像显示延时2秒，然后释放图形(4)为图形窗口加标题程序代码closeclearallx=0:0.05:1;y=f(x);q=length(x);n=input('请输入伯恩斯坦多项式的次数n=[n1,n2,...](各分量是单增的)n=');m=length(n);fori=1:1:mplot(x,y,'b')%绘制连续函数图像holdonforj=1:1:q%计算伯恩斯坦多项式在xj处的值Bjifx(j)==0B(j)=f(0);elseifx(j)==1B(j)=f(1);elsep=(1-x(j))^n(i);%基函数赋初值B(j)=0;fork=0:1:n(i)B(j)=B(j)+f(k/n(i))*p;%计算n(i)次伯恩斯坦多项式的值p=(n(i)-k)*x(j)*p/((k+1)*(1-x(j)));%基函数计算endendendplot(x,B,'r')%绘制n(i)次伯恩斯坦多项式的图像pause(2);holdoffendtitle('伯恩斯坦多项式逼近连续函数')运行结果请输入伯恩斯坦多项式的次数n=[n1,n2,...](各分量是单增的)n=[1,2,3,4,,5,6]程序运行后所得到的图像如图3-1所示。图3-1伯恩斯坦多项式逼近连续函数的动画函数的最佳平方逼近多项式实验用例设，1.利用MATLAB的符号运算功能，给出准确的一次最佳平方逼近多项式；2.利用MATLAB的数值计算功能，给出近似的一次最佳平方逼近多项式；实验步骤及说明1.在MATLAB命令空间以下命令，可以得到法方程的系数矩阵G以及右端向量d。symsx%设x为符号变量G(1,1)=int(x^0,0,1)%调用定积分计算函数int，计算格兰姆矩阵的元素G(1,2)=int(x^1,0,1)G(2,1)=G(1,2)G(2,2)=int(x^2,0,1)d(1)=int(sqrt(1+x^2),0,1)%计算右端向量的元素d(2)=int(x*sqrt(1+x^2),0,1)其结果为G=[1，1/2][1/2，1/3]d=[1/2*2^(1/2)-1/2*log(2^(1/2)-1)，2/3*2^(1/2)-1/3]2.再键入解符号方程组的命令[a1,a2]=solve('a1+a2/2=1/2*2^(1/2)-1/2*log(2^(1/2)-1)','a1/2+a2/3=2/3*2^(1/2)-1/3','a1','a2')可得到a1=-2*2^(1/2)+2-2*log(2^(1/2)-1)，a2=5*2^(1/2)-4+3*log(2^(1/2)-1)这样就得到了一次最佳平方逼近多项式3.调用MATLAB关于定积分数值计算函数quad，同样可以得到数值化的法方程系数矩阵G以及右端向量d，其命令如下：G(1,1)=quad(‘x.^0’,0,1)G(1,2)=quad(‘x.^1’,0,1)G(2,1)=G(1,2)G(2,2)=quad(‘x.^2’,0,1)d(1)=quad(‘sqrt(1+x.^2)’,0,1)d(2)=quad(‘x.*sqrt(1+x.^2)’,0,1)其结果为G=d=1.00000.50001.14780.50000.33330.60954.再键入解数值方程组的命令a=G\d’可得到a=0.93430.4269这样也得到了近似的一次最佳平方逼近多项式希尔伯特矩阵的病态性实验用例对于，权函数，函数系，计算在上的最佳平方逼近多项式时，需要求解法方程。这里，为阶的希尔伯特矩阵而维列向量d的第个分量。通过数值实验，说明法方程的病态性。即，当右端列向量的某个分量出现微小扰动时，其法方程的解会出现较大的差异。实验步骤及说明1.在MATLAB命令空间键入命令G=Hilb(4)，生成一个4阶的希尔伯特矩阵；键入命令d=[1,1,1,1]’，生成一个4维列向量。2.键入命令x=G\d，给出法方程的“精确解”x。3.对d的第1个分量d(1)施加扰动，解法方程，即键入以下命令d(1)=1-0.01x1=G\d4.对d的第4个分量d(4)施加扰动，并再解法方程，即键入以下命令d(4)=1+0.01x2=G\d5.将三个解向量拼成一个矩阵，进行比较，即键入命令[x,x1,x2]其结果为ans=-4.0000-4.1600-5.560060.000061.200078.0000-180.0000-182.4000-224.4000140.0000141.4000169.4000从这一结果不难发现，当第1，第4个分量分别施加0.01的扰动之后，其法方程的解在第3个分量上，竟然相差44.4。这表明，希尔伯特矩阵对右端列向量的取值相当“敏感”。多项式拟合模型的选取实验用例对于数据表xi12345yi2.345.85.46.8分别用一次，二次，三次，四次多项式来拟合这些数据点，并通过作图，找出哪一种拟合多项式对这些数据点的拟合效果最好。实验步骤及说明1.将数据点分别存入向量x和y，即在MATLAB命令空间键入命令x=[1,2,3,4,5];y=[2.3,4,5.8,5.4,6.8];2.利用MATLAB多项式拟合函数polyfit以及多项式求值函数polyval，求出一次多项式拟合函数的系数，并计算一次拟合多项式在x处的函数值。可在MATLAB命令空间键入以下命令p1=polyfit(x,y,1);y1=polyval(p1,x);3.将图形窗口分为4个，利用subplot命令，将数据散点图与一次拟合多项式图像画在第一个窗口。subplot(221)plot(x,y,’*’)holdonplot(x,y1)4.与前两步方法类似地，可完成数据散点图与二次、三次、四次拟合多项式图像。p2=polyfit(x,y,2);y2=polyval(p2,x);subplot(222)plot(x,y,’*’)holdonplot(x,y2)p3=polyfit(x,y,3);y3=polyval(p3,x);subplot(223)plot(x,y,’*’)holdonplot(x,y3)p4=polyfit(x,y,4);y4=polyval(p4,x);subplot(224)plot(x,y,’*’)holdonplot(x,y4)5.这些命令执行后，可得到图3-2。很明显，四次拟合多项式对给定数据点的拟合效果更好一些。图3-2各次拟合多项式的拟合效果比较非线性拟合模型的选取实验用例在某化学反应里，根据实验所得生成物浓度与时间的关系如下表时间(t分)12345678浓度(y10-3)4.006.408.008.809.229.509.709.86时间(t分)910111213141516浓度(y10-3)10.0010.2010.3210.4210.5010.5510.5810.60求浓度y与时间t的拟合曲线y=F(t)。实验步骤及说明1.将数据点分别存入向量x和y，即在MATLAB命令空间键入命令t=1:16;y=［4.00，6.40，8.00，8.80，9.22，9.50，9.70，9.86，10.00，10.20，10.32，10.42，10.50，10.55，10.58，10.60］;2.作数据点的散点图。在MATLAB命令空间键入命令plot(t,0.001*y,’o’)holdon3.图像表明：开始时浓度增加较快，后来逐渐减弱，到一定时间就基本稳定在一个水平。即当t→∞时y趋于某个数，故有一水平渐近线；另外当t=0时，反应末开始，浓度为0。根据这些特点，可设想该数据的数学模型y=F(t)为双曲线型为了确定待定参数，只需对数据(1/ti，1/yi)作一次多项式拟合。键入命令a=polyfit(1./t,1./(0.001*y),1)可得到a=162.722580.1745于是，拟合函数为。4.作双曲拟合曲线的图像。在MATLAB命令空间键入u=1:0.5:16;v=1./(a(1)./u+a(2));plot(u,v,’r’)得到图3-3。从图像上可以看出，拟合效果一般。5.从图像上分析，还可以用指数模型来进行拟合，即拟合曲线为这里：a1<0，当t→∞时，y→a2；当t→0时，y→0，且t增加时y增加，与给出数据规律相同。对该拟合曲线函数取对数，有为了确定待定参数，只需对数据(1/ti，lnyi)作一次多项式拟合。键入命令a=polyfit(1./t,log(0.001*y),1)可得到a=-1.0567-4.4807再键入a(2)=exp(a(2)),得到a=-1.05670.0113于是，拟合函数为。6.再作数据点与指数拟合曲线的图像。在MATLAB命令空间键入figureplot(t,0.001*y,’o’)holdonu=1:0.5:16;v=a(2)*exp(a(1)./u);plot(u,v,’r’)得到图3-4。从图像上可以看出，该模型的拟合效果更好一些。图3-3数据点与双曲拟合图像图3-4数据点与指数拟合图像实验题1.设有下述实验数据表xi1.001.251.501.752.00yi5.105.796.537.458.46(1)分析该实验数据，选择拟合曲线的数学模型；(2)给出拟合函数；(3)作散点图和拟合曲线图，观察拟合效果。2.求函数在指定区间上的最佳一次平方逼近多项式：(1)；(2)3.已知实验数据如下：xi1925313844yi19.032.349.0373.397.8(1)用最小二乘法求形如的经验公式；(2)作数据散点图和拟合曲线图，观察拟合效果。4.下表给出了美国1910年到2000年的人口数，用插值与拟合的方法分别预报2010年美国的人口数。t(年)19101920193019401950y(百万)91.972105.711123.203131.669150.697t(年)19601970198019902000y(百万)179.323203.212226.505249.633281.4225.若点在插值节点之间，用插值函数来预测处的函数近似值的方法，称之为内插法；否则称之为外推法。人口的预测问题如果用插值法如果用插值或拟合的方法，则是外推方法。一般来说，外推法通常是不可靠的。通过下述实例来证实这一观点。设，构造在[2，4]上的拉格朗日插值多项式，对于较大的n，在比[2，4]大的区间上同时画及的图像。数值积分牛顿—科茨公式的数值稳定性实验用例当时，牛顿—科茨系数中的出现了负值，在这种情况下，用牛顿—科茨公式计算定积分的近似值，通常在数值上是不稳定的。1.用牛顿—科茨公式求的近似值；2.根据计算结果，说明在的情况下，其计算结果不可靠。编程思想1.为被积函数创建函数文件f.m，其代码如下：functiony=f(x)y=1/(1+x^2);2.用8行9列矩阵cotes存放牛顿—科茨系数表，用0表示表中的空格。3.为使程序具有一般性，用实型变量a，b分别表示积分区间的左右端点，这里，a=-4，b=4。用n表示牛顿—科茨公式的阶数，这里n=1,2,…,8。4.牛顿—科茨公式的计算算法如下：（1）对于n=1,2,…,8，做以下操作①s←0，h←(b-a)/n②对于k=1,2,…,n+1，做以下操作xk←a+(k-1)h，yk←f(xk)，s←s+cotes(n,k)yk③s←s(b-a)④输出n,s程序代码%创建牛顿—科茨系数矩阵,其中0表示空格cotes=[1/2,1/2,0,0,0,0,0,0,0;1/6,2/3,1/6,0,0,0,0,0,0;1/8,3/8,3/8,1/8,0,0,0,0,0;7/90,16/45,2/15,16/45,7/90,0,0,0,0;19/288,25/96,25/144,25/144,25/96,19/288,0,0,0;41/840,9/35,9/280,34/105,9/280,9/35,41/840,0,0;751/17280,3577/17280,1323/17280,2989/17280,2989/17280,1323/17280,3577/17280,751/17280,0;989/28350,5888/28350,-928/28350,10496/28350,-4540/28350,10496/28350,-928/28350,5888/28350,989/28350];a=input('请输入积分下限a=');b=input('请输入积分上限b=');disp('阶数牛顿—科茨值')forn=1:8s=0;%牛顿—科茨求积公式的初值h=(b-a)/n;%步长fork=1:n+1x(k)=a+(k-1)*h;%节点y(k)=f(x(k));%节点处的函数值s=s+cotes(n,k)*y(k);%计算牛顿—科茨求积公式的和ends=s*(b-a);%牛顿—科茨求积公式的值fprintf('n=%1dI=%5.4f\n',n,s);%输出阶数与牛顿—科茨值end运行结果请输入积分下限a=-4请输入积分上限b=4阶数牛顿—科茨值n=1I=0.4706n=2I=5.4902n=3I=2.2776n=4I=2.2776n=5I=2.3722n=6I=3.3288n=7I=2.7997n=8I=1.9411由，在MATLAB命令空间键入命令2*atan(4)，可得到其值大约为2.6516。从上述程序输出的近似值不难发现，当n>5之后，其计算效果越来越差。这表明，高阶的牛顿—科茨公式由于其数值稳定性差而不宜采用。复化求积思想的动画演示实验用例如果，按照定积分定义，有，（）这实际上是一个最简单的复化求积公式1.用数值实验说明(不妨取）；2.设计一幅动画，摸拟的过程。编程思想1.为被积函数创建函数文件f.m，其代码如下：functiony=f(x)y=exp(x);2.程序的主要算法如下：（1）输入积分下限a，积分上限b以及积分区间的等分数组n（2）m←length(n)（3）对于i=1,2,…,m，做以下操作①作被积函数图像，固化图形②si←0，h←(b-a)/ni③对于j=1,2,…,ni，做以下操作(ⅰ)uj←a+(j-1)h，vj←f(uj)，uj+1←a+jh，si←si+vjh(ⅱ)用黄色填充第j个矩形④图像显示延时2秒，释放图形(4)输出积分区间等分数与阶梯形面积值程序代码closeclearalla=input('请输入积分区间左端点a=');b=input('请输入积分区间右端点b=');n=input('请输入对积分区间等分的正整数数组n=[n1,n2,...](各分量是单增的)n=');m=length(n);%绘制阶梯形逼近曲边梯形的动画,计算阶梯形面积x=a:0.0005:b;y=f(x);fori=1:mplot(x,y,'b')%作被积函数图像holdons(i)=0;h=(b-a)/n(i);forj=1:n(i)u(j)=a+(j-1)*h;v(j)=f(u(j));u(j+1)=a+j*h;s(i)=s(i)+v(j)*h;%计算阶梯形面积fill([u(j),u(j),u(j+1),u(j+1)],[0,v(j),v(j),0],'y')%用黄色填充矩形区域endpause(2)%让图像显示延时2秒holdoffendtitle('阶梯形逼近曲边梯形')%为图像加上标题table=[n',s'];disp('等分数阶梯形面积')disp(table)%输出积分区间等分数与阶梯形面积值运行结果请输入积分区间左端点a=0请输入积分区间右端点b=2请输入对积分区间等分的正整数数组n=[n1,n2,...](各分量是单增的)n=[2,4,8,16,32]等分数阶梯形面积2.00003.71834.00004.92438.00005.623716.00005.998132.00006.1915程序运行后所得到的图像如图4-1所示。图4-1阶梯形面积逼近曲边梯形面积自动选取步长的复化梯形求积算法实验用例对于给定的误差限，利用复化梯形求积的递推公式以及它的截断误差事后估计式编程计算定积分的近似值。编程思想1.为被积函数创建函数文件f.m，其代码如下：functiony=f(x)ifx==0y=1;elsey=sin(x)/x;end2.自动选步长的复化梯形求积算法如下：(1)输入a，b，eps(2)h←b-a，T1←h(f(a)+f(b))/2(3)反复做以下操作①u←h/2，H←0，x←a+u②当x<b时，反复做以下操作(ⅰ)H←H+f(x),x←x+h③T2←(T1+hH)/2④若|T2-T1|<eps，I←T2+(T2-T1)/3，break⑤h←u，T1←T2(4)输出I程序代码a=input('请输入积分下限a=');b=input('请输入积分上限b=');eps=input('请输入误差限eps=');h=b-a;%初始步长T1=h*(f(a)+f(b))/2;%初始梯形值while1u=h/2;%步长两分H=0;%函数值累和变量初始化x=a+u;%初始求积节点whilex<bH=H+f(x);%节点处函数值累加x=x+h;%初始求积endT2=(T1+h*H)/2;%梯形递推ifabs(T2-T1)<epsI=T2+(T2-T1)/3;%带补偿的梯形值breakendh=u;%准备下一轮梯形递推的初值T1=T2;endfprintf('定积分近似值I=%10.8f\n',I)%输出定积分近似值运行结果请输入积分下限a=0请输入积分上限b=1请输入误差限eps=10^(-8)定积分近似值I=0.94608307龙贝格数表实验用例对于，构造其龙贝格下三角数表。编程思想1.为被积函数创建函数文件f.m，其代码如下：functiony=f(x)y=x^(3/2);2.输入下限a，上限b，积分区间的两分次数k。3.将两分次数存入龙贝格数表romberg的第1列。4.给出积分区间两分0次数（即等分1次）的梯形值。5.利用复化梯形公式计算ti(i=2,…,k+1)的值，ti表示积分区间被两分i-1次（即等分2i-1次）的梯形值。6.将梯形值存入romberg矩阵的第2列。7.利用加速公式生成龙贝格数表中的其余各列（其算法与差商计算类似）。8.显示下三角状的龙贝格数表。程序代码clearall%清除命令空间中所有变量formatlong%设置数据类型为长实型a=input('输入定积分下限a=');b=input('输入定积分上限b=');k=input('输入对积分区间两分次数k=');%构造龙贝格表的第1列,即两分次数half=0:1:k;romberg=[half'];%给出积分区间两分0次，即等分1次的梯形值h=b-a;T=h*(f(a)+f(b))/2;t(1)=T;%利用复化梯形公式计算t(2),t(3),...,t(k+1)fori=2:k+1u=h/2;H=0;x=a+u;whilex<bH=H+f(x);x=x+h;endt(i)=(T+h*H)/2;T=t(i);h=u;end%产生龙贝格表的第2列romberg=[romberg,t'];%利用加速公式，生成龙贝格数表的其余各列fori=2:k+1forj=k+1:-1:1ifj>=it(j)=(4^(i-1)*t(j)-t(j-1))/(4^(i-1)-1);elset(j)=0;endendromberg=[romberg,t