张扬2010212420数值计算方法实验讲解

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

2、<math.h>double Lagra nge2(double *x, double *y, double in put)double output;int i;for (i=0;i<5;i+)if (xi <= in put && xi+1 >= in put)output=yi +(yi+1-yi)*(i nput-xi)/(xi+1-xi);break;return output;2. 分段二次插值double Lagra nge3(double *x,double *y,double u)int i,k=0;double v;for(i=

3、0;i<6;i+)if(u<x1)k=0;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-x k+1)/(xk+2-xk)*(xk+2-xk+1);if(xi<u&&u<=xi+1)&&(fabs(u-xi)v=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-

4、xk)*(xk+1-xk+2)+yk+2*(u-xk)*(u-x k+1)/(xk+2-xk)*(xk+2-xk+1);if (xi<u&&u<=xi+1)&&fabs(u-xi)>fabs(u-xi+1)k=i;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-x k+1)/(xk+2-xk)*(xk+2-xk+1);if(u>x4)k=3;v=yk*(u-xk+1)*(u-xk+2)

5、/(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-x k+1)/(xk+2-xk)*(xk+2-xk+1);return v;void mai n()double x6 = 0.0, 0.1, 0.195, 0.3, 0.401,0.5, y6=0.39894,0.39695,0.39142,0.38138,0.36812,0.35206;double u;sca nf("%lf",&u);拉格朗日插值：double Lagra nge1(double *x

6、, 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;double Lagra nge2(double *x, double *y, double in put) double output;int i;for (i=0;i<5;i+)if (xi <= in put && xi+1 >= in

7、 put)output=yi +(yi+1-yi)*(i nput-xi)/(xi+1-xi);break;return output;运行结果:Q Figure 1File Edit View Insert Tools Desktop Window Help'd u a 电q g it b四、四种插值的比较：从以上对比函数图象可以看出，分段线性插值其总体光滑程度不够。在数学上， 光滑程度的定量描述是函数（曲线）的k阶导数存在且连续，则称该曲线具有 k 阶光滑性。一般情况下，阶数越高光滑程度越好。分段线性插值具有零阶光滑性， 也就是不光滑。二次插值就是较低次数的多项式而达到较高阶光滑性

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

9、中，a0，a1,a2,an为待求未知数，n为多项式的最高次幕，由此，该问题 化为求的极值问题。由多元函数求极值的必得到:.j=0,1,cJ锻 羽2 口口#0 i=0-z)a7 =°In mm&址严九二血*=D 7tOi,n这是一个关于a0，a1,a2,an的线性方程组，用矩阵表示如下:mE £rrtNtJ=Omj-0mi=0m.J=or=0求出未知数矩阵(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,'*')xl

10、abel 'x轴'ylabel 'y轴'title ' 散点图hold onm=6;n=3;A=zeros( n+1); for j=1: n+1 for i=1: n+1 for k=1:m+1 A(j,i)=A(j,i)+x(k)A(j+i-2) endendend;B=0 0 0 0;for j=1: n+1for i=1:m+1B(j)=B(j)+y(i)*x(i)A(j-1)endendB=B'a=i nv (A)*B;x=-1.0:0.0001:2.0;z=a(1)+a (2)*x+a (3)*x42+a *x.A3;plot(x,z

11、)legend('离散点','y=a(1)+a(2)*x+a(3)*x.A2+a*x.A3') title(' 拟合图')、运行的结果:Q Figure 1五、经验总结:在运筹学、统计学、逼近论和控制论中，最小二乘法都是很重要的求解方法。例 如，它是统计学中估计回归参数的最基本方法。 在实际问题中，怎样由测量的数 据设计和确定“最贴近”的拟合曲线？关键在选择适当的拟合曲线类型，有时根据专业知识和工作经验即可确定拟合曲线类型；在对拟合曲线一无所知的情况 下，不妨先绘制数据的粗略图形，或许从中观测出拟合曲线的类型；更一般地， 对数据进行多种曲线类型的

12、拟合，并计算均方误差，用数学实验的方法找出在最 小二乘法意义下的误差最小的拟合函数。实验三矩阵的特征值和特征向量、算法设计思想:设A是n阶方阵，若存在数，和n维非零向量x使关系式A XX（1）成立，则称这样的数，为方阵A的特征值；非零向量X称为A的对应于 特征值'的特征向量.将（1）式改写成（A ?uE） x = 0（ 2）得到一个含n个未知数n个方程的齐次线性方程组，它有非零解的充分 必要条件是其系数行列式一A XE =0,、即扎 a2ana21 a22 九a2n=0 A J < Aan1an2ann _九这个以，为未知数的一元 n次方程称为方阵 A的特征方程，其左端 A- &

13、#39;E是，的n次多项式，记作fC）,称为方阵A的特征多项式，显然, A的特征值就是特征方程的根；在复数范围内，n阶方阵有n个特征值。 从定义1不难看出，若非零列向量xi是方阵A的对应于特征值i的特征 向量，则kxi （ k=0为常数）也是A的对应于i的特征向量，即对应于 一个特征值有无穷多个特征向量。反之，不同的特征值所对应的特征向 量不相等，即一个特征向量只能属于一个特征值。由上面的讨论得到方阵A的特征值与特征向量的求法：1）写出方阵A的特征方程|A-入日=0，它的根就是A的全部特征值。2）对每个特征值'0，齐次线性方程组（2）的每一个非零解都是A的对 应于0的特征向量，只要求出

14、（2）的一个基础解系，它们的线性组合（为非零向量时）就是 A的对应于0的全部的特征向量。程序清单:clc;clear;close;A=4,1,-1;16,-2,-2;2,-3,-1;X,B=eig(A)、运行结果:A =3-1-220-22-1-1>> V,D=eig(A)0.7276 -0.57740.62300.4851 -0.5774 -0.24170.4851 -0.57740.74391.0000 000 0.000000 0 1.0000四、所得图形kmax(y(k)(k)(Si (k)、x = y /max(x )x(k+1) = Ay(k)01(1, 1, 1)(1

15、0, 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.475122)66.054548(1,0.716216, -0.24363

16、9)(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, 4.287265, -1.499172)116.00167(1,0.714345,

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

18、个分点对应的函数值。龙贝格求积的基本思路与计算步骤：用复化梯形公式，取区间长度为h，复化梯形公式的值 T(h)与原积分值|ff(x)dx之间存在关系：lf =T(h) a/2 a2h4 |该式称为欧拉-麦克劳林求和公式。上式还表明，用T(h)近似lf的截断误差为: RT(h)=亦2 a2h4O(h2)龙贝格求积算法计算步骤如下：T(0) T° 出Tm 。=臂f(a)+f(b)(1)计算 T° :2T T(2)对分区间a，b并计算(1)1(0) b - aIoIo 2 2ij) . b _ a2iI。和 I。：r a b T (i)1 T(fE T。=2To其中' 为

19、新分点的函数值之和。(0)(i)计算T1与T1 :22T T0(0)2 I 0- I 0(i)T12I122 -1利用外推公式：2k (i '1) i 2kJ 1 k 詔2k (i 1)i2 T0 _ T022k -1Tk22k -1k=1,2,川,n)，(i =0,1,2,Hl,m-k)，直至求(5) l f TT (°)_ T ( mm -1肚名是否真，其中耳为给定的精度。若真，则判断(0)m，否则重复、的计算。(i )k排成三角数表To(O)ToT1(0)To(2)T1T2(0)ToT1(2)T2T3(0)数表沿竖向和斜向都是收敛于If，即当m趋于无穷，T絆与T(T均收

20、敛于lf二、程序清单：1 .第一题：用梯形公式和辛浦生公式计算积分edx，并估计误差。(计算取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

21、 << " the value is " << i << en dl;cout << " Do ne ! " << endl;return 0 ;2 1第二题：用龙贝格法求积分I =-：e=dx，要求误差不超过10。(计算取6位小数)fun ctio n R,quad,err,h=romber(f,a,b, n,tol)M=1;h=b-a;err=1;J=0;R=zeros(4,4);R(1,1)=h*(feval('f,a)+feval('f,b)/2;while(err>

22、;tol)&(J<n )|(J<4)J=J+1;h=h/2;s=0;for p=1:Mx=a+h*(2*p-1);s=s+feval('f,x);endR(J+1,1)=R(J,1)/2+h*s;M=2*M;for K=1:JR(J+1,K+1)=R(J+1,K)+(R(J+1,K)-R(J,K)/(4AK-1);enderr=abs(R(J,J)-R(J+1,K+1);% 取绝对值。endquad=R(J+1,J+1);三、运行结果：测试数据：1 =J亡h必10-5用上述程序计算积分"的值，使其误差不超过解：先用m文件先定义一个名为f.m的函数文件。fu

23、n cti on y=f(x);y=2/sqrt(pi)*exp(x);%exp(x) 以 e 为底的指数。建立一个主程序lbg.mclcclearromberf 0, 1,7, 10八(-6)然后在MATLAB令窗口运行上述主程序，即：>> lbg计算结果如下。ans =2.097800 001.97911.939500 01.94901.93891.93890 01.94141.93891.93891.938901.93951.93891.93891.93891.9389呂 fL1.9389由此可知：'''四、所得图形：(1)复化梯形积分1 functi

24、on Lt J = agui_trapz (. fnamej d2fnamee )2 %fnameA被积函数，为函数fnane的二阶导函熟aj分别为上=界，軌精度3 -y= abs (f eval (d2fname, a: 1 e-51 b);4 - jn=max (y);5 -h=abs (sqrt (12*e/ (b-a). /m);6 -n=ceil (b-a) A)7 -h= (b-a)/n;S-fa=f evala);9 -ft-f eval (fnaniSj b)；10 -f-feval (fnamej a+h: h: b-h+0. 001 *h);11 - t=h* (0. 5*

25、 (fa+fb) +sum (f);121312 ICommand WindowFile Edit Depug Desklop WJndow HelpWarning: MATLAB Tooltok Path Cache is out of date and is not being usedType J help t o olb oz_p at h_ c ache * for more info.To get startedj select HATLAB 丑ejp or Dmmg from the Help menuThe elejuent type "nanerf nust be

26、terminated by the Hatching end-tw</name>*.Could not parse the file: c:lab7toolboxccslirikccslinkinfo.xml>> for mart long>> t=agui_traps(inline( m.*exp(x) inline( (x+2).*exp(x)'),1,2,5e-8)n =7019t =E38905CL2F23022>> t=agui_traps (inline (? 4. / Cl+x* "2厂),inline (? (-

27、8+24. *i. *2), / (l-Ht. "2), 31 )j Oj lj 5e-8)n =3652t =»l(2)复化辛普生求积公式123456789101112131415161718function s = aguisimpsonCfname, d4fnajne, a, b, e )%fname为被积函数，d4fname为函数fname的四阶导函数，a, b分别为上下界，e为精度- y=abs(feval(d4fname, a:le-5:b);- m=jnax (y);- h=abs(2880*e/(b-a)./m). A (0. 25);- n=ceil (b

28、-a) /h)- h=(b-a) /n;一 fa=f eval (fname, a);- fb=feval(fname, b);- s=fa-fb;- x=a;- for i=l:n一x=x+h/2;s=3+4eval(fnamc3 x);- x=x+h/2;s=s+2*feval(fname, x);-end- s=s*h/6;(3)龙贝格积分1234 -5 -6 -789 10 11 - 1213 -14 -15 -16 -17 -18function r = agui_rbg( fname, a$b )%fname为被积函数，a, b分别为下上界e=le-6;i=l；j=l；h=b-a;

29、T (i, l)=h/2*(feval (fname, a)+feval (fname b);T(i+1, l)=T(i, 1 )/2+sum(feval (fname, a+h/2:h:b-h/2+0. 001*h)*h/2;Ki+1, j+l)=4Aj*T(i+l, j)/(4Aj-l)-T(i, j)/(4Aj-l);while abs(T(i+1, i+l)"T(i, i)>ei=i+l;h=h/2;T (i+1, 1)=T (i, l)/2+sum(feval (fname, a+h/2: h: b-h/2+0. 001*h)*h/2;for j=l:iT(i+1, j+l)=4*j*T(i+l, j)/(4*j-l)-T(i, j)/(4*j-l);endendTr=T(i+l, j+1);end-J