数值分析实验报告资料

机电工程学院机械工程陈星星6720150109《数值分析》课程设计实验报告实验一函数插值方法一、问题提出对于给定的一元函数)(xfy的n+1个节点值(),0,1,,jjyfxjn。试用Lagrange公式求其插值多项式或分段二次Lagrange插值多项式。数据如下：（1）jx0.40.550.650.800.951.05jy0.410750.578150.696750.901.001.25382求五次Lagrange多项式5L()x，计算(0.596)f,(0.99)f的值。（提示：结果为(0.596)0.625732f,(0.99)1.05423f）实验步骤：第一步：先在matlab中定义lagran的M文件为拉格朗日函数代码为：function[c,l]=lagran(x,y)w=length(x);n=w-1;l=zeros(w,w);fork=1:n+1v=1;forj=1:n+1if(k~=j)v=conv(v,poly(x(j)))/(x(k)-x(j));endendl(k,:)=v;endc=y*l;end第二步：然后在matlab命令窗口输入：x=[0.40.550.650.80,0.951.05];y=[0.410750.578150.696750.901.001.25382];p=lagran(x,y)回车得到：P=121.6264-422.7503572.5667-377.2549121.9718-15.0845由此得出所求拉格朗日多项式为p（x）=121.6264x5-422.7503x4+572.5667x3-377.2549x2+121.9718x-15.0845第三步：在编辑窗口输入如下命令：x=[0.40.550.650.80,0.951.05];y=121.6264*x.^5-422.7503*x.^4+572.5667*x.^3-377.2549*x.^2+121.9718*x-15.0845;plot(x,y)命令执行后得到如下图所示图形，然后x=0.596;y=121.6264*x.^5-422.7503*x.^4+572.5667*x.^3-377.2549*x.^2+121.9718*x-15.084y=0.6257得到f（0.596）=0.6257同理得到f（0.99）=1.0542（2）jx1234567jy0.3680.1350.0500.0180.0070.0020.001试构造Lagrange多项式6L()x，和分段三次插值多项式，计算的(1.8)f,(6.15)f值。（提示：结果为(1.8)0.164762f,(6.15)0.001266f）实验步骤：第一步定义function[c,l]=lagran(x,y)w=length(x);n=w-1;l=zeros(w,w);fork=1:n+1v=1;forj=1:n+1if(k~=j)v=conv(v,poly(x(j)))/(x(k)-x(j));endendl(k,:)=v;endc=y*l;end定义完拉格朗日M文件第二步：然后在matlab命令窗口输入：x=[1234567];y=[0.3680.1350.0500.0180.0070.0020.001];p=lagran(x,y)回车得到：由此得出所求拉格朗日多项式为p（x）=0.0001x6-0.0016x5+0.0186x4-0.1175x3+0.4419x2-0.9683x+0.9950第三步：在编辑窗口输入如下命令：x=[1234567];y=0.0001*x.^6-0.0016*x.^5+0.0186*x.^4-0.1175*x.^3+0.4419*x.^2-0.9683*x+0.9950;plot(x,y)命令执行后得到如下图所示图形，然后x=1.8;y=4304240283865561*x.^6/73786976294838206464-7417128346304051*x.^5/4611686018427387904+223*x.^4/12000-2821*x.^3/24000+994976512675275*x.^2/2251799813685248-19367*x./20000+199/200y=0.1648得到f（1.8）=0.1648同理得到f（6.15）=0.0013实验结论：插值是在离散数据的基础上补插连续函数，使得这条连续曲线通过全部给定的离散数据点，它是离散函数逼近的重要方法，利用它可通过函数在有限个点处的取值状况，估算出函数在其他点处的近似值。实验二函数逼近与曲线拟合一、问题提出从随机的数据中找出其规律性，给出其近似表达式的问题，在生产实践和科学实验中大量存在，通常利用数据的最小二乘法求得拟合曲线。在某冶炼过程中，根据统计数据的含碳量与时间关系，试求含碳量与时间t的拟合曲线。t(分)05101520253035404550554(10)y01.272.162.863.443.874.154.374.514.584.024.64第一步先写出线性最小二乘法的M文件functionc=lspoly(x,y,m)n=length(x);b=zeros(1:m+1);f=zeros(n,m+1);fork=1:m+1f(:,k)=x.^(k-1);enda=f'*f;b=f'*y';c=a\b;c=flipud(c);第二步在命令窗口输入：lspoly([0,5,10,15,20,25,30,35,40,45,50,55],[0,1.27,2.16,2.86,3.44,3.87,4.15,4.37,4.51,4.58,4.02,4.64],2)回车得到：ans=-0.00240.20370.2305即所求的拟合曲线为y=-0.0024x2+0.2037x+0.2305在编辑窗口输入如下命令：x=[0,5,10,15,20,25,30,35,40,45,50,55];y=-0.0024*x.^2+0.2037*x+0.2305;plot(x,y)命令执行得到如下图实验结论：分析复杂实验数据时,常采用分段曲线拟合方法。利用此方法在段内可以实现最佳逼近,但在段边界上却可能不满足连续性和可导性。分段函数的光滑算法,给出了相应的误差分析.给出了该方法在分段曲线拟合中的应用方法以及凸轮实验数据自动分段拟合。实验四线方程组的直接解法一、问题提出给出下列几个不同类型的线性方程组，请用适当算法计算其解。1、设线性方程组123456789104231210000865365010042213210310215131194426167332386857172635021342530116101191734212246271392012400183248631xxxxxxxxxx5123234613381921(1,1,0,1,2,0,3,1,1,2)Tx2、设对称正定阵系数阵线方程组12345678424024000221213206411418356200216143323218122410394334411142202531011421500633421945xxxxxxxx(1,1,0,2,1,1,0,2)Tx3、三对角形线性方程组123456789104100000000141000000001410000000014100000000141000000001410000000014100000000141000000001410000000014xxxxxxxxxx7513261214455*(2,1,3,0,1,2,3,0,1,1)Tx实验步骤：列主元高斯消去法的matlab的M文件程序function[x,det,index]=Gauss(A,b)%求线形方程组的列主元Gauss消去法，其中，%A为方程组的系数矩阵；%b为方程组的右端项；%x为方程组的解；%det为系数矩阵A的行列式的值；%index为指标变量，index=0表示计算失败，index=1表示计算成功。[n,m]=size(A);nb=length(b);%当方程组行与列的维数不相等时，停止计算，并输出出错信息。ifn~=merror('TherowsandcolumnsofmatrixAmustbeequal!');return;end%当方程组与右端项的维数不匹配时，停止计算，并输出出错信息ifm~=nberror('ThecolumnsofAmustbeequalthelengthofb!');return;end%开始计算，先赋初值index=1;det=1;x=zeros(n,1);fork=1:n-1%选主元a_max=0;fori=k:nifabs(A(i,k))a_maxa_max=abs(A(i,k));r=i;endendifa_max1e-10index=0;return;end%交换两行ifrkforj=k:nz=A(k,j);A(k,j)=A(r,j);A(r,j)=z;endz=b(k);b(k)=b(r);b(r)=z;det=-det;end%消元过程fori=k+1:nm=A(i,k)/A(k,k);forj=k+1:nA(i,j)=A(i,j)-m*A(k,j);endb(i)=b(i)-m*b(k);enddet=det*A(k,k);enddet=det*A(n,n);%回代过程ifabs(A(n,n))1e-10index=0;return;endfork=n:-1:1forj=k+1:nb(k)=b(k)-A(k,j)*x(j);endx(k)=b(k)/A(k,k);end然后在命令窗口输入A=[42-3-1210000;86-5-3650100;42-2-132-1031;0-215-13-1194;-426-167-3323;86-8571726-35;02-13-425301;1610-11-917342-122;462-713920124;00-18-3-24-863-1];b=[51232346133819-21];gauss(A,b)ans=1.0000-1.00000.00001.00002.00000.00003.00001.0000-1.00002.0000(1,1,0,1,2,0,3,1,1,2)Tx高斯-约当消去法maltab的M文件程序function[x,flag]=Gau_Jor(A,b)%求线形方程组的列主元Gauss-约当法消去法，其中，%A为方程组的系数矩阵；%b为方程组的右端项；%x为方程组的解；[n,m]=size(A);nb=length(b);%当方程组行与列的维数不相等时，停止计算，并输出出错信息。ifn~=merror('TherowsandcolumnsofmatrixAmustbeequal!');return;end%当方程组与右端项的维数不匹配时，停止计算，并输出出错信息ifm~=nberror('ThecolumnsofAmustbeequ