2023年数值计算B大作业.doc

课 程 设 计课程名称： 数值计算B 设计题目： 数值计算B大作业 学 号： 姓 名： 完毕时间： 题目一：多项式插值某气象观测站在8：00（AM）开始每隔10分钟对天气作如下观测，用三次多项式插值函数（Newton）逼近如下曲线，插值节点数据如上表，并求出9点30分该地区的温度（x=10）x12345678y22.523.324.421.7025.228.524.825.4二、数学原理假设有n+1个不同的节点及函数在节点上的值（x，y），……（x，y），插值多项式有如下形式： （1）其中系数（i=0,1,2……n）为特定系数，可由插值样条（i=0,1,2……n）拟定根据均差的定义，把x当作[a,b]上的一点,可得f(x)= f（）+f[]（）f[x, ]= f[]+f[x,] （）……f[x, ，…x]= f[x, ，…x]+ f[x, ，…x](x-x)综合以上式子，把后一式代入前一式，可得到： f(x)= f[]+f[]（）+ f[]（）（）+…+ f[x, ，…x]（）…(x-x)+ f[x, ，…x,]= N（x）+其中N（x）= f[]+f[]（）+ f[]（）（）+…+ f[x, ，…x]（）…(x-x) （2）= f(x)- N（x）= f[x, ，…x,] （3） =（）…(x-x)Newton插值的系数（i=0,1,2……n）可以用差商表达。

一般有[] （k=0,1,2，……，n ） （4）把（4）代入（1）得到满足插值条件N（i=0,1,2，……n）的n次Newton插值多项式N（x）=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 = 21.7313A = 22.5000 0 0 0 24.4000 0.9500 0 0 25.2023 0.4000 -0.1375 0 24.8000 -0.2023 -0.1500 -0.0021C = -0.0021 -0.1187 1.4521 21.1688 L = - x^3/480 - (19*x^2)/160 + (697*x)/480 + 3387/160四、 结果分析和讨论 对于不超过三次的插值多项式，x假如选取1,3,5,7这三个点可以得到较好的三次插值多项式L=-0.0021x^3-0.1187x^2+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 0 1.27 2.16 2.86 3.44 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 * -0.0024 0.2037 0.2305>> plot(x,y,'r')四、 结果分析和讨论 通过最小二乘法得到的曲线拟合多项式是p=(-0.0024x^2+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,且r≠n. 求解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=l\b %输出向量yx=u\y %输出向量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 -。