课 程 设 计课程名称： 数值计算B 设计题目： 数值计算B大作业 学 号： S316030014 姓 名： 周帅宇 完成时间： 2021年11月3日 题目一：多项式插值某气象观测站在8：00AM开始每隔10分钟对天气作如下观测，用三次多项式插值函数Newton逼近如下曲线，插值节点数据如上表，并求出9点30分该地区的温度x=10。x12345678y二、数学原理假设有n+1个不同的节点及函数在节点上的值x，y，x，y，插值多项式有如下形式： 1其中系数i=0,1,2n为特定系数，可由插值样条i=0,1,2n确定。根据均差的定义，把x看成a,b上的一点,可得f(x)= f+ffx, = f+fx, fx, ，x= fx, ，x+ fx, ，x(x-x)综合以上式子，把后一式代入前一式，可得到： f(x)= f+f+ f+ fx, ，x(x-x)+ fx, ，x,= Nx+其中Nx= f+f+ f+ fx, ，x(x-x) 2= f(x)- Nx= fx, ，x, 3 =(x-x)Newton插值的系数i=0,1,2n可以用差商表示。一般有 k=0,1,2，n 4把4代入1得到满足插值条件Ni=0,1,2，n的n次Newton插值多项式Nx=f+f+f+f.其中插值余项为： 介于之间。三、程序设计function y,A,C,L=newdscg(X,Y,x,M) % y为对应x的值，A为差商表，C为多项式系数，L为多项式% X为给定节点，Y为节点值，x为待求节点n=length(X); m=length(x); % n为X的长度for t=1:m z=x(t); A=zeros(n,n);A(:,1)=Y;s=0.0; p=1.0; q1=1.0; c1=1.0; for j=2:n for i=j:n A(i,j)=(A(i,j-1)- A(i-1,j-1)/(X(i)-X(i-j+1); end q1=abs(q1*(z-X(j-1);c1=c1*j; end C=A(n,n);q1=abs(q1*(z-X(n);for k=(n-1):-1:1C=conv(C,poly(X(k); d=length(C);C(d)=C(d)+A(k,k);end y(k)= polyval(C, z); %输出y值endL(k,:)=poly2sym(C); %输出多项式 syms M,X=1,3,5,7;Y=22.5,24.4,25.2,24.8;x=10; y,A,C,L=newdscg(X,Y,x,M)y =A = 22.5000 0 0 0 24.4000 0.9500 0 0 25.2000 0.4000 -0.1375 0C = L = - x3/480 - (19*x2)/160 + (697*x)/480 + 3387/160四、 结果分析和讨论 对于不超过三次的插值多项式，x如果选取1,3,5,7这三个点能够得到较好的三次插值多项式L=-0.0021x3-0.1187x2+1.4521x+21.1688。当x=10时，也即9点30分时的温度为21.7317度，结果分析知此值应是偏小的。对于选取不同的插值节点，能够得到不同的插值多项式，误差也不尽相同。五、 完成题目的体会与收获 牛顿插值法的重要一点就是对插值节点的选取，通过此题的编程很好的加深了对其概念的理解以及提高了应用牛顿插值法的能力，学会了运用Matlab软件对牛顿插值法相关问题进行编程求解，对Matlab计算方法与程序编辑更加熟悉。使我对这类问题的理解有了很大的提升。题目二：曲线拟合在某钢铁厂冶炼过程中，根据统计数据的含碳量与时间关系，试用最小二乘法拟合含碳量与时间t的拟合曲线，并绘制曲线拟合图形。t(分)0 5 10 15 20 25 30 35 40 45 50 55 3.87 4.15 4.37 4.51 4.58 4.02 4.64 二、 数学原理从整体上考虑近似函数同所给数据点(i=0,1,m)误差(i=0,1,m)的大小，常用的方法有以下三种：一是误差(i=0,1,m)绝对值的最大值，即误差 向量的范数；二是误差绝对值的和，即误差向量r的1范数；三是误差平方和的算术平方根，即误差向量r的2范数；前两种方法简单、自然，但不便于微分运算 ，后一种方法相当于考虑 2范数的平方，因此在曲线拟合中常采用误差平方和来 度量误差(i=0，1，m)的整体大小。数据拟合的具体作法是：对给定数据 (i=0,1,，m)，在取定的函数类中,求,使误差(i=0,1,m)的平方和最小，即 =从几何意义上讲，就是寻求与给定点(i=0,1,m)的距离平方和为最小的曲线。函数称为拟合 函数或最小二乘解，求拟合函数的方法称为曲线拟合的最小二乘法。在曲线拟合中，函数类可有不同的选取方法。三、程序设计 x=0,5,10,15,20,25,30,35,40,45,50,55; y=0,1.27,2.16,2.86,3.44,3.87,4.15,4.37,4.51,4.58,4.02,4.64*10(-4); p=polyfit(x,y,2)p = 1.0e-04 * plot(x,y,r)四、 结果分析和讨论 通过最小二乘法得到的曲线拟合多项式是p=(-0.0024x2+0.2037x+0.2305)*10-4。利用Matlab软件调用最小二乘法函数即可得到多项式拟合函数，由于数据较少得到的拟合曲线不够光滑，可能会存在一定的误差。拟合曲线总体呈现增函数趋势，与数据较为吻合。五、 完成题目的体会与收获 曲线拟合较常用的方法之一就包括最小二乘法，因此能够熟练使用Matlab进行数据的曲线拟合变得尤为重要。通过完成此题的拟合过程，对于最小二乘法曲线拟合的操作更加的熟练，能够较好的完成各类数据的拟合。题目三：线性方程组求解分别利用LU分解；平方根法与改良平方根法；追赶法求解以下几个不同类型的线性方程组，并与准确值比拟结果，分析误差产生原因。1设线性方程组 二、数学原理 将A分解为一个下三角矩阵L和一个上三角矩阵U，即：A=LU,其中 L=， U= 解Ax=b 的问题就等价于要求解两个三角形方程组： Ly=b,求y; Ux=y,求x。 设A为非奇异矩阵,且有分解式A=LU, L为单位下三角阵,U为上三角阵。 L,U的元素可以有n步直接计算定出。用直接三角分解法解Ax=b要求A的所有顺序主子式都不为零的计算公式： , ,i=2,3,，n. 计算U的第r行，L的第r列元素(i=2,3,n)： ， i=r,r+1,n; , i=r+1,n,且rn. 求解Ly=b，Ux=y的计算公式； 三、程序设计function f1=LUresolve(a,b);n,n=size(a); % 行数z=size(b) % b的行数l=; % l矩阵u=; % u矩阵for j=1:1:n u(1,j)=a(1,j);endfor i=2:1:n l(i,1)=a(i,1)/u(1,1);endfor i=2:1:(n-1) sum1=0; for k=1:1:(i-1) sum1=l(i,k)*u(k,i)+sum1; end u(i,i)=a(i,i)-sum1; sum2=0; sum3=0; for j=(i+1):1:n for k=1:1:(i-1) sum2=sum2+l(i,k)*u(k,j); sum3=sum3+l(j,k)*u(k,i); end u(i,j)=a(i,j)-sum2; l(j,i)=(a(j,i)-sum3)/u(i,i); endendfor i=1:1:(n-1) l(i,i)=1; l(i,n)=0;endl(n,n)=1;sum4=0;for k=1:1:(n-1) sum4=sum4+l(n,k)*u(k,n);endu(n,n)=a(n,n)-sum4;l %输出l矩阵u %输出u矩阵y=lb %输出向量yx=uy %输出向量x a=4 2 -3 -1 2 1 0 0 0 0 8 6 -5 -3 6 5 0 1 0 0 4 2 -2 -1 3 2 -1 0 3 1 0 -2 1 5 -1 3 -1 1 9 4 -4 2 6 -1 6 7 -3 3 2 3 8 6 -8 5 7 17 2 6 -3 5 0 2 -1 3 -4 2 5 3 0 1 16 10 -11 -9 17 34 2 -1 2 2 4 6 2 -7 13 9 2 0 12 4 0 0 -1 8 -3 -24 -8 6 3 -1; b=5;12;3;2;3;46;13;38;19;-21; LUresolve(a,b)z =