数值分析第五版计算实习题
数值分析计算实习题二
《数值分析》计算实习题二算法设计方案1．主要计算步骤：计算函数f(x,y)在拟合所需的节点处的函数值。
将各拟合节点（x i，y j）分别带入非线性方程组0.5 cos t + u + v + w – x = 2.67t + 0.5 sin u + v + w – y = 1.070.5t + u + cos v + w – x =3.74t + 0.5u + v + sin w – y =0.79解非线性方程组得解向量（t ij，u ij，v ij，w ij）。
对数表z(t,u)进行分片二次代数插值，求得对应（t ij，u ij）处的值，即为f(x i，y j) 的值。
对上述拟合节点分别进行x，y最高次数为k（k=0,1,2,3…）次的多项式拟合。
每次拟合后验证误差大小，直到满足要求。
2．求解非线性方程组选择Newton迭代法，迭代过程中需要求解线性方程组，选择选主元的Doolittle分解法。
3．对z(t，u)进行插值选择分片二次插值。
4．拟合基函数φr(x)ψs(y)选择为φr(x)=x r，ψs(y)=y s。
拟合系数矩阵c通过连续两次解线性方程组求得。
一．源程序#include "stdio.h"#include "stdlib.h"#include "math.h"void Doolittle(double *A,int n,int *M)//功能说明：对n阶矩阵A进行选主元的Doolittle分解//参数说明：A：欲进行分解的方阵，同时也是返回参数，分解后的结果// 存储于A中// n：方阵A的维数// M；（返回参数）n维向量，记录选主元过程中行交换的次序{int i,j,k,t;double *s;double Maxs,temp;s=(double*) calloc(n,sizeof(double));for(k=0;k<n;k++){for(i=k;i<n;i++){s[i]=A[i*n+k];for(t=0;t<k;t++) s[i]-= A[i*n+t] * A[t*n+k];}Maxs=abs(s[k]); M[k]=k;for(i=k+1;i<n;i++){if(Maxs<abs(s[i])){Maxs=abs(s[i]);M[k]=i;}}if(M[k]!=k){for(t=0;t<n;t++){temp=A[k*n+t];A[k*n+t]=A[M[k]*n+t];A[M[k]*n+t]=temp;}temp=s[k];s[k]=s[M[k]];s[M[k]]=temp;}if(Maxs<(1e-14)){s[k]=1e-14;printf("%.16e方阵奇异\n",Maxs);}A[k*n+k]=s[k];for(j=k+1;(j<n)&&(k<n-1);j++){for(t=0;t<k;t++) A[k*n+j]-=A[k*n+t]*A[t*n+j];A[j*n+k]=s[j]/A[k*n+k];}}}void Solve_LUEquation(double* A,int n,double* b,double* x) //功能说明：解方程LUx＝b，其中L、U共同存储在A中//参数说明：A：经Doolittle分解后的方阵// n：方阵A的维数// b：方程组的右端向量// x：（返回参数）方程组的解向量{int i,t;for(i=0;i<n;i++){x[i]=b[i];for(t=0;t<i;t++) x[i]-=A[i*n+t]*x[t];}for(i=n-1;i>-1;i--){for(t=i+1;t<n;t++) x[i]-=A[i*n+t]*x[t];x[i]/=A[i*n+i];}}void Transpose(double *A,int m,int n,double* AT)//功能说明：求m×n阶矩阵A的转置AT//参数说明：A：已知m×n阶矩阵// m：A的行数// n：A的列数// AT：（返回参数）A的转置矩阵（n×m）{int i,j;for(i=0;i<m;i++)for(j=0;j<n;j++) AT[j*m+i]=A[i*n+j];}void Solve_LEquation(double* A,int n,double* B,double* x,int m) //功能说明：解线性方程组Ax＝B，该函数可对系数矩阵相同// 而右端向量不同的多个方程组同时求解。
最新数值分析(第五版)计算实习题第五章作业
数值分析第五章第一题：LU分解法：建立m文件function h1=zhijieLU(A,b)%h1各阶主子式的行列式值[n n]=size(A);RA=rank(A);if RA~=ndisp('请注意：因为A的n阶行列式h1等于零，所以A不能进行LU分解。
A的秩RA如下：')RA,h1=det(A);returnendif RA==nfor p=1:nh(p)=det(A(1:p,1:p));endh1=h(1:n);for i=1:nif h(1,i)==0disp('请注意：因为A的r阶主子式等于零，所以A不能进行LU分解。
A的秩RA和各阶顺序主子式h1依次如下：')h1;RAreturnendendif h(1,i)~=0disp('请注意：因为A的r阶主子式都不等于零，所以A能进行LU分解。
A的秩RA和各阶顺序主子式h1依次如下：')for j=1:nU(1,j)=A(1,j);endfor k=2:nfor i=2:nfor j=2:nL(1,1)=1;L(i,i)=1;if i>jL(1,1)=1;L(2,1)=A(2,1)/U(1,1);L(i,1)=A(i,1)/U(1,1);L(i,k)=(A(i,k)-L(i,1:k-1)*U(1:k-1,k))/U(k,k);elseU(k,j)=A(k,j)-L(k,1:k-1)*U(1:k-1,j);endendendendh1;RA,U,L,X=inv(U)*inv(L)*bendend输入：>> A=[10 -7 0 1;-3 2.099999 6 2;5 -1 5 -1;2 1 0 2];>> b=[8;5.900001;5;1];>> h1=zhijieLU(A,b)输出：请注意：因为A的r阶主子式都不等于零，所以A能进行LU分解。
A的秩RA和各阶顺序主子式h1依次如下：RA =4U =10.0000 -7.0000 0 1.00000 2.1000 6.0000 2.30000 0 -2.1429 -4.23810 -0.0000 0 12.7333L =1.0000 0 0 0-0.3000 1.0000 0 00.5000 1.1905 1.0000 -0.00000.2000 1.1429 3.2000 1.0000X =-0.2749-1.32981.29691.4398h1 =10.0000 -0.0000 -150.0001 -762.0001列主元高斯消去法：建立m文件function [RA,RB,n,X]=liezhu(A,b)B=[A b];n=length(b);RA=rank(A);RB=rank(B);zhicha=RB-RA;if zhicha>0disp('请注意：因为RA~=RB，所以方程组无解')returnwarning off MATLAB:return_outside_of_loopendif RA==RBif RA==ndisp('请注意：因为RA=RB，所以方程组有唯一解')X=zeros(n,1);C=zeros(1,n+1);for p=1:n-1[Y,j]=max(abs(B(p:n,p)));C=B(p,:);B(p,:)=B(j+p-1,:);B(j+p-1,:)=C;for k=p+1:nm=B(k,p)/B(p,p);B(k,p:n+1)=B(k,p:n+1)-m*B(p,p:n+1);endendb=B(1:n,n+1);A=B(1:n,1:n);X(n)=b(n)/A(n,n);for q=n-1:-1:1X(q)=(b(q)-sum(A(q,q+1:n)*X(q+1:n)))/A(q,q);endelsedisp('请注意：因为RA=RB<n，所以方程组有无穷多解') endend输入：>> A=[10 -7 0 1;-3 2.099999 6 2;5 -1 5 -1;2 1 0 2];>> b=[8;5.900001;5;1];>> [RA,RB,n,X]=liezhu(A,b),H=det(A)输出：请注意：因为RA=RB，所以方程组有唯一解RA =4RB =4n =4X =0.0000-1.00001.00001.0000H =-762.0001第二题：建立列主元高斯消去法m文件（题一中已有）（1）输入：>> format compact>> A=[3.01 6.03 1.99;1.27 4.16 -1.23;0.987 -4.81 9.34];>> b=[1;1;1];>> [RA,RB,n,X]=liezhu(A,b),h=det(A),C=cond(A)输出：请注意：因为RA=RB，所以方程组有唯一解RA =3RB =n =3X =1.0e+03 *1.5926-0.6319-0.4936h =-0.0305C =3.0697e+04(2)输入：>> A=[3.00 6.03 1.99;1.27 4.16 -1.23;0.990 -4.81 9.34]; >> b=[1;1;1];>> [RA,RB,n,X]=liezhu(A,b),h=det(A)输出：请注意：因为RA=RB，所以方程组有唯一解RA =3RB =3n =3X =119.5273-47.1426-36.8403h =-0.4070第三题：输入：>> clear>> A=[10 7 8 7;7 5 6 5;8 6 10 9;7 5 9 10];>> b=[32 23 33 31]’;>> dA=det(A),lamda=eig(A),Ac2=cond(A,2)输出：dA =1.0000lamda =0.01020.84313.858130.28872.9841e+03下面分析误差性态：建立m文件：function Acp=pjwc(A,jA,b,jb,p)%Acp矩阵A的p条件数cond%pjwc：p范数解的误差性态分析%jA是A的近似矩阵jA=A+δA，jb=b+δbAcp=cond(A,p);dA=det(A);X=A\b;deltaA=jA-A;pndA=norm(deltaA,p);deltab=jb-b;pndb=norm(deltab,p);if pndb>0jX=A\jb;Pnb=norm(b,p);pnjx=norm(jX,p);deltaX=jX-X;pnjdX=norm(deltaX,p);jxX=pnjdX/pnjX;pnX=norm(X,p);xX=pnjdX/pnX;pndb=norm(deltab,p);xAb=pndb/pnb;pnbj=norm(jb,p);xAbj=pndb/pnbj;Xgxx=Acp*xAb;endif pndA>0jX=jA\b;deltaX=jX-X;pnX=norm(X,p);pnjdX=norm(deltaX,p);pnjX=norm(jX,p);jxX=pnjdX/pnjX;xX=pnjdX/pnX;pnjA=norm(jA,p);pnA=norm(A,p);pndA=norm(deltaA,p);xAbj=pndA/pnjA;xAb=pndA/pnA;Xgxx=Acp*xAb;endif (Acp>50)&(dA<0.1)disp('请注意：AX=b是病态的，A的p条件数Acp，A的行列式值dA，解X，近似解jX，解的相对误差xX，解的相对误差估计Xgxx，b或A的相对误差xAb依次如下：')Acp,dA,X',jX',xX',jxX',Xgxx',xAb',xAbj'elsedisp('请注意：AX=b是良态的，A的p条件数Acp，A的行列式值dA，解X，近似解jX，解的相对误差xX，解的相对误差估计Xgxx，b或A的相对误差xAb依次如下：')Acp,dA,X',jX',xX',jxX',Xgxx',xAb',xAbj'end输入：>> jA=[10 7 8.1 7.2;7.08 5.04 6 5;8 5.98 9.89 9;6.99 5 9 9.98];>> jb=b;p=2;>> Acp=pjwc(A,jA,b,jb,p)输出：请注意：AX=b是良态的，A的p条件数Acp，A的行列式值dA，解X，近似解jX，解的相对误差xX，解的相对误差估计Xgxx，b或A的相对误差xAb依次如下：Acp =2.9841e+03dA =1.0000ans =1.0000 1.0000 1.0000 1.0000 ans =-9.5863 18.3741 -3.2258 3.5240 xX =10.4661jxX =0.9842Xgxx =22.7396xAb =0.0076xAbj =0.0076Acp =2.9841e+03第四题：（1）输入：建立m文件:for n=2:6a=hilb(n);pnH(n-1)=cond(a,inf);endpnHn=2:6;plot(n,pnH);可见条件数随着n的增大而急剧增大（2）输入：>> n=2;H=hilb(n);>> x=(linspace(1,1,n))';>> b=H*x;>> [RA,RB,n,X]=gauss(H,b)输出：请注意：因为RA=RB，所以方程组有唯一解RA =2RB =2n =2X =1.00001.0000输入：>> r=b-H*X,deltax=X-x输出：r =deltax =1.0e-15 *0.4441-0.6661输入：>> n=3;H=hilb(n);>> x=(linspace(1,1,n))';>> b=H*x;>> [RA,RB,n,X]=gauss(H,b)>> r=b-H*X,deltax=X-x输出：X =1.00001.00001.0000r =1.0e-15 *0.2220deltax =1.0e-13 *-0.02000.1221-0.1255>> n=4;H=hilb(n);>> x=(linspace(1,1,n))'; >> b=H*x;>> [RA,RB,n,X]=gauss(H,b) >> r=b-H*X,deltax=X-xX =1.00001.00001.00001.0000r =1.0e-15 *-0.4441-0.1110deltax =1.0e-12 *-0.02220.2485-0.59800.3886>> n=5;H=hilb(n);>> x=(linspace(1,1,n))'; >> b=H*x;>> [RA,RB,n,X]=gauss(H,b) >> r=b-H*X,deltax=X-xX =1.00001.00001.00001.00001.0000r =1.0e-15 *0.22200.1110deltax =1.0e-11 *-0.00350.0524-0.19370.2591-0.1148>> n=6;H=hilb(n);>> x=(linspace(1,1,n))'; >> b=H*x;>> [RA,RB,n,X]=gauss(H,b) >> r=b-H*X,deltax=X-xX =1.00001.00001.00001.00001.0000r =1.0e-15 *0.22200.1110deltax =1.0e-11 *-0.00350.0524-0.19370.2591-0.1148>> n=7;H=hilb(n);>> x=(linspace(1,1,n))'; >> b=H*x;>> [RA,RB,n,X]=gauss(H,b) >> r=b-H*X,deltax=X-xX =1.00001.00001.00001.00001.00001.0000r =1.0e-15 *0.22200.1110deltax =1.0e-09 *-0.00080.0219-0.14820.3854-0.42540.1677>> n=8;H=hilb(n);>> x=(linspace(1,1,n))'; >> b=H*x;>> [RA,RB,n,X]=gauss(H,b) >> r=b-H*X,deltax=X-xX =1.00001.00001.00001.00001.00001.00001.00001.0000r =1.0e-15 *-0.2220-0.1110-0.1110deltax =1.0e-06 *-0.00000.0018-0.02360.1279-0.34420.4870-0.34660.0978>> n=9;H=hilb(n);>> x=(linspace(1,1,n))'; >> b=H*x;>> [RA,RB,n,X]=gauss(H,b) >> r=b-H*X,deltax=X-xX =1.00001.00001.00001.00001.00001.00001.00001.00001.0000r =1.0e-15 *0.44410.2220-0.22200.22200.2220-0.1110deltax =1.0e-04 *-0.00000.0002-0.00280.0197-0.07220.1471-0.16870.1017-0.0251>> n=10;H=hilb(n);>> x=(linspace(1,1,n))';>> b=H*x;>> [RA,RB,n,X]=gauss(H,b)>> r=b-H*X,deltax=X-xX =1.00001.00001.00001.00000.99991.00030.99961.00040.99981.0000r =1.0e-15 *0.4441-0.22200.2220-0.11100.11100.1110deltax =1.0e-03 *-0.00000.0001-0.00230.0205-0.09740.2669-0.43690.4214-0.22090.0485一、填空题(每空1分，共40分)1、中国共产党是_中国工人阶级___的先锋队，同时是中国人民和中华民族的先锋队，是_______中国特色社会主义事业___ 的领导核心。
数值分析课程实验设计——数值积分实习题
数值分析——数值积分实习题管理科学与工程学院 学号：1120140500 姓名：彭洋洋 一、计算实习题1.用不同数值方法计算积分：049xdx =-⎰.（1）取不同的步长h ，分别用复合梯形及复合辛普森求积计算积分，给出误差中关于h 的函数，并与积分精确值比较两个公式的精度，是否存在一个最小的h ，使得精度不能再被改善？ （2）用龙贝格求积计算完成问题（1） （3）用自适应辛普森积分，使其精度达到10-4解答：（1）取不同的步长，采用不同的公式，比较精度过程如下： 1.1 复合梯形公式及复合辛普森公式求解复合梯形公式：11*[()2()()]2n n k k hT f a f x f b -==++∑误差关于h 的函数：2(2)()**()12n a b R f h f ξ-=复合辛普森公式：111/201*[()4()2()()]6n n n k k k k hS f a f x f x f b --+===+++∑∑误差关于h 的函数：4(4)()*(/2)*()180n a bR f h f η-=1.2 复合梯形公式及复合辛普森公式Matlab 程序（2）用龙贝格求积计算完成问题（1） 2.1 龙贝格求积算法龙贝格求积公式也称为逐次分半加速法。
它是在梯形公式、辛普森公式和柯特斯公式之间的关系的基础上，构造出一种加速计算积分的方法。
作为一种外推算法，它在不增加计算量的前提下提高了误差的精度。
24133n n n S T T =- 21611515n n n C S S =- 26416363n n n R C C =-1221/201()22n n n k k h T T f x -+==+∑ ()(1)()11(4*)/(41)k m k k mm m m T T T +--=-- 1,2,...k = 2.2 龙贝格求积Matlab 程序画图程序设计①得到关于n各种公式求积的图表如下：对于梯形公式、辛普森公式n代表份数，龙贝格公式n表示从1开始的序列号②关于步长h 的各种公式求积的图表如下其中龙贝格序列步长()/2k h b a =-：观察两幅图表h 越小，精度越高。
数值分析第五章实习题答案
数值分析第五章实习题答案数值分析第五章实习题答案数值分析是一门研究如何使用计算机来解决数学问题的学科。
在数值分析的学习过程中，实习题是非常重要的一部分，通过实习题的练习，可以帮助我们巩固所学的知识，并且提高我们的解题能力。
本文将为大家提供数值分析第五章实习题的答案，希望对大家的学习有所帮助。
第一题：求下列方程的一个正根，并用二分法和牛顿法分别计算根的近似值。
方程：x^3 - 3x + 1 = 0解答：首先，我们可以通过绘制函数图像来初步估计方程的根的范围。
根据图像，我们可以大致确定根在区间[0, 2]之间。
接下来，我们使用二分法来计算根的近似值。
根据二分法的原理，我们将区间[0, 2]等分为两部分，然后判断根在哪一部分。
不断重复这个过程，直到找到根的近似值。
具体计算过程如下：- 将区间[0, 2]等分为两部分，得到中点x = 1。
- 计算方程在x = 1处的函数值f(1) = -1。
- 根据函数值的正负性，我们可以确定根在区间[1, 2]之间。
- 将区间[1, 2]等分为两部分，得到中点x = 1.5。
- 计算方程在x = 1.5处的函数值f(1.5) = 1.375。
- 根据函数值的正负性，我们可以确定根在区间[1, 1.5]之间。
- 重复以上步骤，直到找到根的近似值。
最终得到根的近似值为x ≈ 1.365。
接下来，我们使用牛顿法来计算根的近似值。
牛顿法是一种迭代法，通过不断逼近根的位置来计算根的近似值。
具体计算过程如下：- 选择初始近似值x0 = 1。
- 计算方程在x = 1处的函数值f(1) = -1。
- 计算方程在x = 1处的导数值f'(1) = 4。
- 利用牛顿法的迭代公式x1 = x0 - f(x0)/f'(x0)，我们可以得到x1 ≈ 1.333。
- 重复以上步骤，直到找到根的近似值。
最终得到根的近似值为x ≈ 1.365。
通过二分法和牛顿法，我们分别得到了方程x^3 - 3x + 1 = 0的一个正根的近似值为x ≈ 1.365。
数值分析计算实习作业一
数值分析计算实习题一学号：：院系：2015年11月5日一、分析1.1算法分析题目要求求出：1）特征值从小到大排列的最小特征值1λ和最大特征值501λ。
2）特征值中模最小的特征值s λ。
3）靠近一组数k μ的一组特征值k i λ。
4）矩阵A 的条件数cond(A)2。
5）行列式detA 。
解决方法：1）若将所有行列式按模的大小排列则模最大的特征值一定是1λ和501λ中的一个，因此利用幂法求出模最大的特征值1m λ。
然后利用带原点平移的幂法，将系数矩阵变为1m A I λ-即将所有特征值都减去1m λ，则特征值按大小顺序排列的次序不变，模最大的特征值依然在整个排列的两端，再用一次幂法得到模最大的特征值21=m m λλλ-，其中λ为带原点平移的幂法求出的特征值，最后两个特征值1m λ、2m λ比较大小，大的为501λ，小的为1λ。
2）因为s λ为按模最小的特征值，因此用反幂法可求的其特征值。
3）因为k i λ靠近数k μ，因此k i k λμ-一定是所有的k λμ-中模最小的，因此可利用带原点平移的反幂法求出特征值k i λ，此时的系数矩阵变为k A I μ-。
4）条件数cond(A)2为模最小的特征值与模最大的特征值的比的绝对值，因此利用1和2中求出的1m λ和s λ可解出条件数。
5）可对矩阵A 进行LU 分解，即A LU =则det()det()det()A L U =⨯，又因为矩阵L 对角线元素为1，则det()L =1，所以det()det()A U =，U 为上三角阵，行列式为对角线元素的乘积，因此可得A 的行列式。
1.2程序分析1.2.1 因为A 为拟三角阵，储存时零元素不储存，因此将矩阵A 压缩为5*501的矩阵CA 的带元素ij a =C 中的元素1,i j s j c -++ 程序中A[5][501]即为压缩后的矩阵。
1.1.2 程序中的B[5][501]为过渡矩阵，在幂法迭代、反幂法迭代以及LU 分解中均用矩阵B 来计算，计算之间对B 进行适当的赋值。
数值分析计算实习题
插值法1.下列数据点的插值x 0 1 4 9 16 25 36 49 64y 0 1 2 3 4 5 6 7 8可以得到平方根函数的近似，在区间[0,64] 上作图.(1)用这9 个点作8 次多项式插值Ls(x).(2)用三次样条( 第一边界条件)程序求S(x).从得到结果看在[0,64] 上，哪个插值更精确；在区间[0,1] 上，两种插值哪个更精确解：(1) 拉格朗日插值多项式，求解程序如下syms x l;x1=[0 1 4 9 16 25 36 49 64]; y1=[0 1 2 3 4 5 6 7 8]; n=length(x1);Ls=sym(0);for i=1:nl=sym(y1(i));for k=1:i-1l=l*(x-x1(k))/(x1(i)-x1(k));endfor k=i+1:n l=l*(x-x1(k))/(x1(i)-x1(k));endLs=Ls+l; endLs=simplify(Ls) % 为所求插值多项式Ls(x).输出结果为Ls = -/*xA2+95549/72072*x-1/00*xA8-2168879/0*xA4+19/0*xA7+657859/*xA3+33983/ 0*xA5-13003/00*xA6(2) 三次样条插值，程序如下x1=[0 1 4 9 16 25 36 49 64];y1=[0 1 2 3 4 5 6 7 8];x2=[0:1:64];y3=s plin e(x1,y1,x2);p=po Iyfit(x2,y3,3); % 得到三次样条拟合函数S=p(1)+p(2)*x+ p(3)*x^2+p(4)*xA3 % 得到S(x) 输出结果为：S =/6464-2399/88*x+/1984*xA2+2656867/624*xA3⑶ 在区间［0,64］上，分别对这两种插值和标准函数作图,Plot(x2,sqrt(x2),'b',x2,y2,'r',x2,y3,'y')蓝色曲线为y="X函数曲线，红色曲线为拉格朗日插值函数曲线，黄色曲线为三次样条插值曲线可以看到蓝色曲线与黄色曲线几乎重合，因此在区间[0,64] 上三次样条插值更精确。
数值分析计算实习题
插值法1.下列数据点的插值x 0 1 4 9 16 25 36 49 64y 0 1 2 3 4 5 6 7 8可以得到平方根函数的近似，在区间[0,64]上作图.(1)用这9个点作8次多项式插值Ls(x).(2)用三次样条(第一边界条件)程序求S(x).从得到结果看在[0,64]上，哪个插值更精确；在区间[0,1]上，两种插值哪个更精确？解：(1)拉格朗日插值多项式，求解程序如下syms x l;x1=[0 1 4 9 16 25 36 49 64];y1=[0 1 2 3 4 5 6 7 8];n=length(x1);Ls=sym(0);for i=1:nl=sym(y1(i));for k=1:i-1l=l*(x-x1(k))/(x1(i)-x1(k));endfor k=i+1:nl=l*(x-x1(k))/(x1(i)-x1(k));endLs=Ls+l;endLs=simplify(Ls) %为所求插值多项式Ls(x).输出结果为Ls =-24221063/63504000*x^2+95549/72072*x-1/3048192000*x^8-2168879/43545600 0*x^4+19/283046400*x^7+657859/10886400*x^3+33983/152409600*x^5-13003/2395 008000*x^6(2)三次样条插值，程序如下x1=[0 1 4 9 16 25 36 49 64]; y1=[0 1 2 3 4 5 6 7 8]; x2=[0:1:64];y3=spline(x1,y1,x2);p=polyfit(x2,y3,3); %得到三次样条拟合函数 S=p(1)+p(2)*x+p(3)*x^2+p(4)*x^3 %得到S(x) 输出结果为：S =23491/304472833/8*x+76713/*x^2+6867/42624*x^3(3)在区间[0,64]上，分别对这两种插值和标准函数作图，plot(x2,sqrt(x2),'b',x2,y2,'r',x2,y3,'y')蓝色曲线为y=函数曲线,红色曲线为拉格朗日插值函数曲线，黄色曲线为三次样条插值曲线010203040506070-2020406080100可以看到蓝色曲线与黄色曲线几乎重合，因此在区间[0,64]上三次样条插值更精确。
数值分析第五版计算实习题
弟二草插值法3.卜列数据点的插值可以得到平方根函数的近似，在区间064］上作图。
(1〉用这9个点做8次多项式插值Q x)。
(2)用三次样条(第一边界条件)程岸求S(X)。
从得到结果石在［0.64］ 1：・哪个插值更粘确：在区间［0,1］ I：•两种插值哪个更精确?(1) 8次多项式插值：(1)8次多项式插值:首先建立新的M-file:输入如卜代码(此为拉格朗口插值的功能函数)并保存function f=Language(x,y,x0)%求Li知数据点的拉格朗Fl插值多项式%己知数据点的x坐标向量：x%已知数据点的y坐标向量:y%插值的x坐标：x0%求得的拉格朗H插值多项式或在X0处的插值：fsyms t;ifi(lcngth(x)=length(y))n=length(x);elsedisp(*x和y的维数不相等!)；return;end %检错tbr(i=l:n)i=y(i);fbr(j=1:i-l)l=l*(t-x(j))/(x(i)-x(j))；end;for(j=i-M:n)end;for(j=i+l:n) l=l*(t-x(j))/(x(i)-x(j))； end;simplify(f);if(i==n) if|nargin=3)f=subs(C't\xO);else f=collcct(f);f=vpa(f,6);endendend再建立新的M-file：输入：clear;x=[0 1 49 16 25 36 49 64];y=[0：l：8]；%计算拉格朗口基丞数%计算拉格朗ri插值函数%化简%计算插值点的曲数值%将插值多项式展开%将插值多项式的系数化成6位精度的小数f=Uinguage(x,y) 运行得到f=1.32574*1-381410*t A2+.604294e-1 *t A3+.222972e-3 *t A5-.542921 e-5*t A6+.671268e・7T7・.328063e・9T8・.498071 e-2*t A4 这就是8次多项式插值L s(x)= 1.32574怜.381410*t A2+.604294e-1 *t A3+.222972e-3 *t A5-.542921 e-5*t A6+.671268e-7*t A7-.328063e-9*t A8-.498071 e-2*t A4. (2)三次样条插值：建立新的M-filc：输入：clear;x=[0 I 49 1625 36 4964];尸[0:8];t=[0:0.1:64];Y=t.A(0.5);O=Language(x,y)f= 1,32574*t-.381410*t.A2+.604294e-1 *t.A3+.222972e-3*t.A5-.542921 e・5*(. W+.671268e-7*t.A7-.328063e-9*t.A8-.498071 e-2 *t.A4;S=interp l(x,y,t.'spline,);plol(x,y,o；(・YY.lf.'b'」S'g:');grid;运行程序得到如下图：从结果屮很明显可以看出在[0.64].上.三次样条插值更精确，儿乎与原函数帀合。