张扬2010212420数值计算方法实验

1、 实验一 插值法一、 算法设计思想:在插值计算中，取样点的多少往往会影响所得插值函数优化程度。一般情况下，取样点越多所得插值函数越优化，对应的函数值与标准函数值越接近。在MATLAB中实现分段线性插值，分段二次插值，拉格朗日插值的命令为interp1,其调用格式为: Y1=interp1(X,Y,X1，method)函数根据X,Y的值，计算函数在X1处的值。X,Y是两个等长的已知向量，分别描述采样点和样本值，X1是一个向量或标量，描述欲插值点，Y1是一个与X1等长的插值结果。二、 程序清单：1. 分段线性插值#include#include double Lagrange2(double *x

2、, double *y, double input) double output; int i; for (i=0;i5;i+) if (xi = input) output=yi +(yi+1-yi)*(input-xi)/(xi+1-xi); break; return output;2.分段二次插值double Lagrange3(double *x,double *y,double u) int i,k=0; double v; for(i=0;i6;i+) if(ux1) k=0; v=yk*(u-xk+1)*(u-xk+2)/(xk-xk+1)*(xk-xk+2)+yk+1*(u-

3、xk)*(u-xk+2)/(xk+1-xk)*(xk+1-xk+2)+yk+2*(u-xk)*(u-xk+1)/(xk+2-xk)*(xk+2-xk+1); if(xiu&u=xi+1)&(fabs(u-xi)=fabs(u-xi+1) k=i-1; v=yk*(u-xk+1)*(u-xk+2)/(xk-xk+1)*(xk-xk+2)+yk+1*(u-xk)*(u-xk+2)/(xk+1-xk)*(xk+1-xk+2)+yk+2*(u-xk)*(u-xk+1)/(xk+2-xk)*(xk+2-xk+1); if (xiu&ufabs(u-xi+1) k=i; v=yk*(u-xk+1)*(u-

4、xk+2)/(xk-xk+1)*(xk-xk+2)+yk+1*(u-xk)*(u-xk+2)/(xk+1-xk)*(xk+1-xk+2)+yk+2*(u-xk)*(u-xk+1)/(xk+2-xk)*(xk+2-xk+1); if(ux4) k=3; v=yk*(u-xk+1)*(u-xk+2)/(xk-xk+1)*(xk-xk+2)+yk+1*(u-xk)*(u-xk+2)/(xk+1-xk)*(xk+1-xk+2)+yk+2*(u-xk)*(u-xk+1)/(xk+2-xk)*(xk+2-xk+1); return v;void main() double x6 = 0.0, 0.1, 0

5、.195, 0.3, 0.401, 0.5,y6 = 0.39894,0.39695,0.39142,0.38138,0.36812,0.35206; double u; scanf(%lf,&u);拉格朗日插值：double Lagrange1(double *x, double *y, double xx) int i,j; double *a,yy=0.000; a=new double6; for(i=0;i 6;i+) ai=yi; for(j=0;j 6;j+) if(j!=i) ai*=(xx-xj)/(xi-xj); yy+=ai; delete a; return yy;do

6、uble Lagrange2(double *x, double *y, double input) double output; int i; for (i=0;i5;i+) if (xi = input) output=yi +(yi+1-yi)*(input-xi)/(xi+1-xi); break; return output;三、 运行结果：四、 四种插值的比较：从以上对比函数图象可以看出，分段线性插值其总体光滑程度不够。在数学上，光滑程度的定量描述是函数(曲线) 的k阶导数存在且连续，则称该曲线具有k阶光滑性。一般情况下，阶数越高光滑程度越好。分段线性插值具有零阶光滑性，也就是不光

7、滑。二次插值就是较低次数的多项式而达到较高阶光滑性的方法。总体上分段线性插值具有以下特点：五、经验总结：分段线性插值在计算上具有简洁方便的特点。分段线性插值与二次插值函数在每个小区间上相对于原函数都有很强的收敛性，（舍入误差影响不大），数值稳定性好且容易在计算机上编程实现等优点，但分段线性插值在节点处具有不光滑性的缺点(不能保证节点处插值函数的导数连续)，从而不能满足某些工程技术上的要求。实验二 曲线拟合的最小二乘法一、 算法设计思想：对给定数据 (i=0 ,1,2,3,.,m),一共m+1个数据点，取多项式P(x),使二、函数P(x)称为拟合函数或最小二乘解，令似的 使得 其中，a0，a1,

8、a2,an为待求未知数，n为多项式的最高次幂，由此，该问题化为求的极值问题。由多元函数求极值的必要条件： j=0,1,n得到： j=0,1,n这是一个关于a0，a1,a2,an的线性方程组，用矩阵表示如下：因此，只要给出数据点 及其个数m，再给出所要拟合的参数n，则即可求出未知数矩阵（a0，a1,a2,an）二、程序清单： x=-1.0:0.5:2.0;y=-4.447,-0.452,0.551,0.048,-0.447,0.549,4.552;plot(x,y,*)xlabel x轴ylabel y轴title 散点图hold onm=6;n=3;A=zeros(n+1);for j=1:n

9、+1 for i=1:n+1 for k=1:m+1 A(j,i)=A(j,i)+x(k)(j+i-2) end endend;B=0 0 0 0;for j=1:n+1 for i=1:m+1 B(j)=B(j)+y(i)*x(i)(j-1) endendB=B;a=inv(A)*B;x=-1.0:0.0001:2.0;z=a(1)+a(2)*x+a(3)*x.2+a(4)*x.3;plot(x,z) legend(离散点,y=a(1)+a(2)*x+a(3)*x.2+a(4)*x.3)title(拟合图)三、运行的结果：五、经验总结：在运筹学、统计学、逼近论和控制论中，最小二乘法都是很重要

10、的求解方法。例如，它是统计学中估计回归参数的最基本方法。在实际问题中，怎样由测量的数据设计和确定“最贴近”的拟合曲线？关键在选择适当的拟合曲线类型，有时根据专业知识和工作经验即可确定拟合曲线类型；在对拟合曲线一无所知的情况下，不妨先绘制数据的粗略图形，或许从中观测出拟合曲线的类型；更一般地，对数据进行多种曲线类型的拟合，并计算均方误差，用数学实验的方法找出在最小二乘法意义下的误差最小的拟合函数。实验三 矩阵的特征值和特征向量一、算法设计思想：设是阶方阵，若存在数和维非零向量使关系式 （1）成立，则称这样的数为方阵A的特征值；非零向量称为A的对应于特征值的特征向量将（1）式改写成， （2）得到一

11、个含个未知数个方程的齐次线性方程组，它有非零解的充分必要条件是其系数行列式一、即这个以为未知数的一元次方程称为方阵的特征方程，其左端是的次多项式，记作，称为方阵的特征多项式，显然，的特征值就是特征方程的根；在复数范围内，阶方阵有个特征值。从定义1不难看出，若非零列向量是方阵的对应于特征值的特征向量，则（为常数）也是的对应于的特征向量，即对应于一个特征值有无穷多个特征向量。反之，不同的特征值所对应的特征向量不相等，即一个特征向量只能属于一个特征值。由上面的讨论得到方阵的特征值与特征向量的求法：1）写出方阵的特征方程，它的根就是的全部特征值。2）对每个特征值，齐次线性方程组（2）的每一个非零解都是

12、的对应于的特征向量，只要求出（2）的一个基础解系，它们的线性组合（为非零向量时）就是的对应于的全部的特征向量。二、 程序清单：clc;clear;close; A=4,1,-1;16,-2,-2;2,-3,-1; X,B=eig(A)三、运行结果：A = 3 -1 -2 2 0 -2 2 -1 -1 V,D=eig(A)V = 0.7276 -0.5774 0.6230 0.4851 -0.5774 -0.2417 0.4851 -0.5774 0.7439D = 1.0000 0 0 0 0.0000 00 0 1.0000赞同四、所得图形kmax(y(k)x(k) = y(k)/max(x

13、(k)x(k+1) = Ay(k)01(1, 1, 1)(10, 8, 1)110(1, 0.8, 0.1)(7.2, 5.4, -0.8)27.2(1, 0.75, -0.111111)(6.5, 4.75, -1.222222)36.57(1, 0.730769, -0.203704)(6.230766, 4.499997, -1.407408)46.230766(1, 0.722222, -0.225880)(6.111108, 4.388886, -1.1451767)56.111108(1, 0.718182, -0.237561)(6.054548, 4.336336, -1.47

14、5122)66.054548(1, 0.716216, -0.243639)(6.027024, 4.310808, -1.487278)76.027024(1, 0.715247, -0.246768)(6.013458, 4.298211, -1.483536)86.013458(1, 0.714765, -0.248366)(6.00671, 4.291945, -1.496732)96.00671(1, 0.714525, -0.249177)(6.00335, 4.28825, -1.496354)106.00335(1, 0.714405, -0.249586)(6.00167,

15、4.287265, -1.499172)116.00167(1, 0.714345, -0.239792)(6.00083, 4.286485, -1.499584)126.00083(1, 0.714315, -0.249896)五、经验总结：MATLAB提供了eig来计算矩阵的特征值、特征向量信息。如果再结合使用MATLAB的排序函数等资源，可以综合利用来解决问题。实验四 数值积分一、算法设计思想：辛普生公式,亦即为抛物线法求其近似值,其实质是将曲线弧换为二次抛物线弧,再对抛物线积分。其结果为:由于“曲线”y=f(x)在区间a,b上的平均“高度(”值)为:(如图)代入(1)式得:(2)其中

16、:yi(i=0,1,2,32n)为将区间a,b进行2n等份,得2n+1个分点,依次每个分点对应的函数值。龙贝格求积的基本思路与计算步骤：用复化梯形公式，取区间长度为，复化梯形公式的值与原积分值之间存在关系： 该式称为欧拉-麦克劳林求和公式。上式还表明，用近似的截断误差为：龙贝格求积算法计算步骤如下：(1) 计算：(2) 对分区间并计算和：，其中为新分点的函数值之和。(3) 计算与：，(4) 利用外推公式：，()，()，直至求出。(5) 判断是否真，其中为给定的精度。若真，则，否则重复(2)、(3)、(4)的计算。排成三角数表数表沿竖向和斜向都是收敛于，即当趋于无穷，与均收敛于。二、程序清单：第

17、一题：用梯形公式和辛浦生公式计算积分，并估计误差。（计算取5位小数）int main() int a=1 ;int b=2 ;int n=100 ; double h = 0.01 ;double xk = 0 ;double xk1 = 0 ;double xk12 = 0 ;double i = 0 ;for ( int k =0; k=n-1; k+ )xk = a + k*h ;xk12 = xk + h/2 ;xk1 = xk + h ;i = i + f(xk) + 4*f(xk12) + f( xk1 ) ; i = i * h/6 ; cout the value is i endl; cout Done ! tol)&(Jn)|(J lbg计算结果如下。ans = 2.0978 0 0 0 0 1.9791 1.9395 0 0 0 1.9490 1.9389 1.9389 0 0 1.9414 1.9389 1.9389 1.9