多篇数值分析实践报告_matlab

数值分析实验报告拉格朗日插值法拉格朗日插值法基本原理：通过平面上不同两点可以确定一条直线，这就是拉格朗日线性插值问题，对于不在同一条直线的三个点得到的插值多项式则为抛物线。

拉格朗日插值的基多项式（即基函数）为：ni x x ••x •x x l ij j ji i i ,,2,1,0,)(0 =--=∏≠=有了基函数以后就可以直接构造如下多项式：∑==ni i i n x l x f p 0)()(该多项式就是拉格朗日插值法所求得的插值多项式。

拉格朗日插值法算法：1、根据所给点),(i i y x 的坐标依次写出其差值基函数（用for 循环可以轻易解决）ni x x ••x •x x l ij j ji i i ,,2,1,0,)(0 =--=∏≠=2、将差值基函数与其对应的点的函数值相乘得：n i x l x f i i ,,2,1,0),()( =3、将2中各项累加即得差值多项式：∑==ni i i n x l x f p 0)()(拉格朗日插值法程序：function lagrange(A)%A 为一个只有两行的矩阵，第一行为插值点，第二行为插值点对应的函数值 [m,n]=size(A); f=1;p=0;%两个用到的变量 syms xfor i=1:nf=(x-A(1,i))*f;endfor j=1:ng(j)=f/(x-A(1,j));%求插值基函数的分母h(j)=subs(g(j),x,A(1,j));%求插值基函数的分子s(j)=g(j)/h(j)*A(2,j);%插值基函数endfor k=1:ns(j)=collect(s(j));%合并同类项endfor i=1:np=p+s(i);endfprintf('拉格朗日插值法可得多项式:')collect(p) % 可用lagrange([1 3 6 8;4 6 9 12])调试例：有如下表格中有四个插值点及其对应的函数值，用lagrange插值法写出其三次插值多项式：x 1 3 6 8iy 4 6 9 12i解：在matlab命令窗口输入：lagrange([1 3 6 8;4 6 9 12])可得运行结果：拉格朗日插值法可得多项式:ans =x^3/70 - x^2/7 + (97*x)/70 + 96/35牛顿插值法牛顿插值法基本原理：拉格朗日插值多项式的理论在许多方面都有应用，是很不错的一种方法，但就插值问题而言，如果增加一个插值点，原先计算插值的多项式就没有用了，即拉格朗日插值法的继承性很差，牛顿插值法就很好地解决了这个问题，具有很好的继承性。

函数)(x f 的差商定义为：)(][k k x f x f =kk k k k k x x x f x f x x f --=---111][][],[kk k k k k k k k x x x x f x x f x x x f --=------212112],[],[],,[jk k k j k k j k k j k j k x x x x f x x f x x x f ---+-+----=],,[],,[],,,[111在差商的基础上可得牛顿插值多项式：)())(](,,,[))(](,,[)](,[][)(11010102100100----++--+-+=n n n x x x x x x x x x f x x x x x x x f x x x x f x f x N牛顿插值法算法：1、根据差商的定义求出各阶差商。

2、依次写出：1,,1,0),())((110-=---=-n n x x x x x x n n ω。

3、n ω分别与相应的差商相乘，然后累加即可得：)())(](,,,[))(](,,[)](,[][)(11010102100100----++--+-+=n n n x x x x x x x x x f x x x x x x x f x x x x f x f x N牛顿插值法程序：function newton(A)%A 为两行矩阵，第一行为插值点，第二行为对应插值点的函数值 [m,n]=size(A);B=zeros(n-1);k=1;for i=1:n-1for j=i+1:nZ(j)=(A(2,j)-A(2,j-1))/(A(1,j)-A(1,j-k));%求i阶差商B(i,j)=Z(j);%记住各阶差商endfor j=i+1:nA(2,j)=Z(j);endk=k+1;endf=A(2,1);g=1;syms xfor i=1:n-1g=g*(x-A(1,i));f=f+B(i,i+1)*g;%求牛顿插值多项式endfprintf('牛顿插值法可得多项式:')ffprintf('牛顿插值法得到的多项式合并同类项后为:')collect(f)% newton([1 3 6 8;4 6 9 12])例：有如下表格中有四个插值点及其对应的函数值，用牛顿插值法写出其三次插值多项式：x 1 3 6 8iy 4 6 9 12i解：在matlab命令窗口输入：newton([1 3 6 8;4 6 9 12])可得如下运行结果：牛顿插值法可得多项式:f =3+x+1/70*(x-1)*(x-3)*(x-6)牛顿插值法得到的多项式合并同类项后为:ans =96/35+1/70*x^3-1/7*x^2+97/70*x通过比较对同一题目的运算可知：拉格朗日插值法与牛顿插值法所得的插值多项式是一致的！最小二乘法最小二乘法基本原理：最小二乘法是一种数学优化技术，它通过最小化误差的平方和寻找数据的最佳函数匹配。

利用最小二乘法可以简便地求得未知的数据，并使得这些求得的数据与实际数据之间误差的平方和为最小。

最小二乘法还可用于曲线拟合，这无疑就是我们数值分析所研究的最小二乘法很重要的一块的。

设已知某一组观测数据m i x f x i i ,,2,1)),(,( =，要求在某特定函数类)(x φ寻求一个函数)(x ϕ作为该组数据的近似函数，使得二者在i x 上的残差 m i x f x i i i ,,2,1),()( =-=ϕδ，按某种度量标准为最小，这就是拟合问题。

要求残差i δ按某种度量标准为最小，即要求由残差i δ构成的残差向量T m ],,,[10δδδδ =的某种范数δ为最小，现实生活中，我们用得最多的当然就是2范数了，即∑∑==-==mi i i mi ix f x 021202122})]()([{)(ϕδδ为最小.这种要求误差平方和最小的拟合称为曲线拟合的最小二乘法.就是说，最小二乘法提供了一种数学方法，利用这种方法可以对实验数据实现在最小平方误差意义下的最好拟合。

在曲线拟合中，函数类φ可有不同的选取方法，我们数值分析中用的最多的自然是最简单与我们最熟悉的多项式。

最小二乘法算法：1、 构造法方程组：具体做法是针对已知的点的坐标),(i i y x ，先求∑=inih n x,1,0,2 再求∑=ii n i h h n y x ,1,0, 为所要拟合的多项式的次数。

2、解法方程组：根据1中构造的方程组，调用列主元高斯消元法（具体算法见列主元高斯消元法的试验）。

3、根据2中求出的法方程组的解X 构造拟合多项式：∑+=-111)(h i i x i X 。

最小二乘法程序：function leastsqu(A,n)%A 代表离散点的坐标矩阵，第一行为自变量值，第二行为相对应的函数值 %n 为要拟合的多项式的次数 [a,b]=size(A); B=zeros(1,10); C=zeros(1,10); P=0;Q=0;D=zeros(n+1);k=1;g=0;if b>n^2s=b;elses=n^2;end %一系列用到的变量或者矩阵for i=1:sfor j=1:bP=P+A(1,j)^(i-1);Q=Q+A(1,j)^(i-1)*A(2,j);endB(i)=P;C(i)=Q;P=0;Q=0;end%for循环嵌套求法方程组的各项系数for i=1:n+1for j=i:n+iD(i,k)=B(j);k=k+1;endD(i,n+2)=C(i);k=1;end%for循环嵌套求法方程组增广矩阵f=qiujie(D);syms xdisp('最小二乘拟合多项式为：')for i=1:n+1g=g+f(i)*x^(i-1);endgfunction f=qiujie(B)[m,n]=size(B);for i=2:mfor k=i:mif(abs(B(k,i-1))>abs(B(i-1,i-1)))temp=B(i-1,:);B(i-1,:)=B(k,:);B(k,:)=temp;%换行选主元endl(i,k)=B(k,i-1)/B(i-1,i-1);for j=1:nB(k,j)=B(k,j)-l(i,k)*B(i-1,j);%跌代 endendendy=0;x(m)=B(m,n)/B(m,n-1);for i=m-1:-1:1for j=1:m-iy=y+x(m-j+1)*B(i,n-j);endx(i)=(B(i,n)-y)/B(i,n-m-1+i);%回代计算结果 y=0;endf=x;例：已知数据表如下，试用二次多项式来拟合。

x0 1 2 3 4 5 6iy15 14 14 14 14 15 16i解：在matlab命令窗口输入：A=[0 1 2 3 4 5 6;15 14 14 14 14 15 16];n=2;leastsqu(A,n)得运行结果：最小二乘拟合多项式为：g =(5*x^2)/28 - (25*x)/28 + 209/14列主元高斯消元法列主元高斯消元法基本原理：列主元高斯消元法，是在高斯消元法的基础上改进的一个线性代数中的算法，用来求解线性方程组。

其原理是：设有方程组：⎪⎪⎩⎪⎪⎨⎧=+++=+++=+++n n nn n n n n n n b x a x a x a bx a x a x a b x a x a x a 22112222212111212111首先写出方程组的增广矩阵，然后利用初等行变换来把增广矩阵转换成行阶梯阵。

特别注意的是每步消元时，选出主元所在列绝对值最大的数，并将其换到主元的位置上，然后进行消元计算，这是列主元消元在高斯消元基础上的改进，是保证算法稳定性之所在。

接着将之转换为行阶梯阵，判断解的情况，但为了简便，我们在这就仅考虑方程组中方程个数与变量个数相等的情况，即只有唯一解的方程组的情况。

列主元高斯消元法算法：1、 令1=k 。

2、 选主元：选出k 列中位于对角线及其以下元素绝对值中最大者，并将该最大数所在行与第k 行互换。