计算方法与软件.doc

实验目的作为实践性非常强的课程，安排上机实验的目的，不仅是为了验证教材和授课内容，更重要的是，要通过实验深入理解方法的设计原理与处理问题的技巧，培养自行处理常规数值计算问题的能力和综合运用知识分析、解决问题的能力。1、通过上机实验加深课堂内容的理解。数值分析的主要任务就是研究适合于在计算机上使用的数值计算方法及与此相关的理论。通过编程上机，就可以加深对方法运行过程的理解，同时在编程中领会和理解数值计算方法的计算要领和步骤，体会问题的条件和限制范围，理解一般问题和特殊问题的区别。2、学会对数值计算结果的分析和处理。数值分析实验不只是编写程序得到一个数值结果，我们应在掌握数值计算计算方法的基本原理和思想的同时，注意方法处理的技巧及其与计算机的密切结合，重视误差分析、收敛性及稳定性的讨论。此外，还要注意算法能否在计算机上实现，应避免因数值方法选用不当、程序设计不合理而导致超过计算机的存储能力，或导致计算结果精度不高等。3、要能灵活掌握各种数值计算方法。由于针对同一个问题可以选用不同的数值计算方法，我们要注意各种方法的使用条件。通过上机，比较各种方法间的异同及优缺点，以便更好的使用不同的方法来解决实际问题，使计算机成为我们最好的工具。实验基本要求一、上机前的准备工作1、复习和掌握与本次实验有关的教学内容。2、根据本次实验要求，在纸上编写算法及上机的程序，并经过人工模拟运行检验，减少不必要的错误，提高上机效率。切忌不编程序、不作人工检查就进行程序输入，这只能使上机调试的难度增加，甚至可能带来学习自信心的下降，影响后续课程的学习。二、上机实验步骤1、启动开发环境；2、建立源程序文件，输入源程序；3、编译产生目标程序，连接生成可执行程序，运行程序，输出结果；4、对数值计算结果进行误差分析，讨论数值算法的收敛性与稳定性；5、整理实验报告。三、实验报告实验报告是记录实验工作全过程的技术文档，实验报告的撰写是科学技术工作的一个组成部分。数值分析实验报告包括下列要求：1、 实验原理；2、 实验内容和要求；3、 数值算法描述，包括数据输入、数据处理和数据输出；4、算法的实现（1） 给出具体的计算实例，（2） 经调试正确的源程序清单，（3） 对具体的数值例子给出数值结果；5、计算结果的误差分析，算法的收敛性与稳定性的讨论；6、实验心得。实验项目实验一、误差分析误差问题是数值分析的基础，又是数值分析中一个困难的课题。在实际计算中，如果选用了不同的算法，由于舍入误差的影响，将会得到截然不同的结果。因此，选取算法时注重分析舍入误差的影响，在实际计算中是十分重要的。同时，由于在数值求解过程中用有限的过程代替无限的过程会产生截断误差，因此算法的好坏会影响到数值结果的精度。一、实验目的1、 通过上机编程，复习巩固以前所学程序设计语言及上机操作指令；2、 通过上机计算，了解误差、绝对误差、误差界、相对误差界的有关概念；3、 通过上机计算，了解舍入误差所引起的数值不稳定性。二、算法实例例1.1 用差商求在处导数的近似值。取，=0.000 000 000 000 001和=0.000 000 000 000 000 1分别用MATLAB软件计算，取十五位数字计算。解： 在MATLAB工作窗口输入下面程序a=3;h=0.1;y=log(a+h)-log(a);yx=y/h运行后得yx = 0.32789822822991.将此程序中改为0.000 1，运行后得yx = 0.33332777790385.后者比前者好。再取h = 0.000 000 000 000 001，运行后得yx = 0.44408920985006，不如前者好。取h = 0.000 000 000 000 000 1，运行后得yx = 0，算出的结果反而毫无价值。例1.2 分别求方程组在下列情况时的解，其中.（1）; （2）.解： （1） 首先将方程组化为同解方程，然后在MATLAB工作窗口输入程序 b=2,2;A=1,1;1,1.01; X=Ab运行后输出当时，的解为；（2）同理可得，当时，的解为.例1.3 计算的近似值。解：泰勒级数e ,取，得. （1.1）这是一个无限过程，计算机无法求到精确值。只能在（1.1）取有限项时计算，再估计误差。如果取有限项作为的值必然会有误差，根据泰勒余项定理可知其截断误差为e.如果取（1.1）的前九项，输入程序 n=8;s=1;S =1;fork=1:ns=s*k;S=S+1/s,ends, S,R=3/(s*(n+1)或S1=1+1+1/2+1/(1*2*3)+1/(1*2*3*4)+1/(1*2*3*4*5)+1/(1*2*3*4*5*6)+1/(1*2*3*4*5*6*7)+1/(1*2*3*4*5*6*7*8),R1=3/(1*2*3*4*5*6*7*8*9)运行后结果S =8.267195767195768e-006 R =2.71827876984127因为截断误差为所以e的近似值e2.718 28.例1.4 取作为的四舍五入近似值时，求其绝对误差和相对误差。解：在MATLAB工作窗口输入程序juewu=exp(1)-2.71828运行后输出结果为juewu = 1.828 459 045 505 326e-006例1.5 计算d 的近似值，并确定其绝对误差和相对误差。解 因为被积函数的原函数不是初等函数，故用泰勒级数求之。 , （1.5）这是一个无限过程，计算机无法求到精确值。可用（1.5）的前四项代替被积函数，得d)d=.根据泰勒余项定理和交错级数收敛性的判别定理，得到绝对误差= WU，在MATLAB命令窗口输入计算程序如下：syms xf=1-x2/(1*2*3)+x4/(1*2*3*4*5)-x6/(1*2*3*4*5*6*7)y=int(f,x,0,pi/2),y1=double(y)y11=pi/2-(pi/2)3/(3*3*2)+(pi/2)5/(5*5*4*3*2)-(pi/2)7/(7*7*6*5*4*3*2)inf=int(sin(x)/x,x,0,pi/2) ,infd=double(inf)WU =(pi/2)9/(9*9*8*7*6*5*4*3*2), R =infd-y11因为运行后输出结果为： 1.370 762 168 154 49，=1.370 744 664 189 38，1.750 396 510 491 47e-005， WU= 1.782 679 830 970 664e-005. 所以，的绝对误差为，故d。的相对误差为 k,juecha,xiangcha,xk= liti112(2,2.5,100) 运行后输出计算结果列入表11和表 1-2中。 将算法2的MATLAB调用函数程序的函数分别用y1=15-2*x2和y1=x-(2*x2+x-15)/(4*x+1)代替，得到算法1和算法3的调用函数程序，将其保存，运行后将三种算法的前8个迭代值列在一起（见表 1-1），进行比较.将三种算法的对应的绝对误差和相对误差的值列在一起（见表 1-2），进行比较。表 1-1 例1.6中三种算法的计算结果算 法迭代次数算法1的迭代结果算法2的迭代结果算法3的迭代结果022.000 000 002.000 000 00173.000 000 002.555 555 562-832.142 857 142.500 550 063-13 7632.837 837 842.500 000 064-378 840 3232.246 963 562.500 000 005-2.870 42.246 963 562.500 000 006-1.647 82.321 774 842.500 000 007-5.430 72.657 901 652.500 000 0099-Inf2.500 000 012.500 000 00表 1-2 例1.6中三种算法计算结果的误差算法迭代 次数算法1的误差算法2的误差算法3的误差绝对误差相对误差绝对误差相对误差绝对误差相对误差00.500 000 000.250 000 000.500 000 000.250 000 000.500 000 000.250 000 0014.500 000 000.642 857 140.500 000 000.166 666 670.055 555 600.021 739 13285.500 000 001.030 120 480.357 142 860.1666 666 700.000 550 100.000 219 97313 765.500 000.000 100 020.337 837 840.119 047 620.000 000 060.000 000 024378 840 3261.000 000 010.253 036 440.112 612 620.000 000 000.000 000 0052.870 399 8110.230 287 040.084 345 480061.647 839 0110.178 225 160.076 762 470075.430 746 8010.157 901 650.059 408 390099InfNaN0.000 000 010.000 000 0000 例1.7 求数的近似值。解 （1）直接用MATLAB命令 x=(715)*(sqrt(1+8(-19)-1)运行后输出结果x = 0.问题出现在两个相近的数与相减时，计算机运行程序sqrt(1+8(-19)-1运行后输出结果 ans = 0.由于计算机硬件只支持有限位机器数的运算，因此在计算中可能引入和传播舍入误差.因为有效数字的严重损失，导致输出的结果为0，计算机不能再与数继续进行真实的计算，所以，最后输出的结果与的精确值不符。（2）如果化为，再用MATLAB命令 x=(715)*( (8(-19)/(sqrt(1+8(-19)+1)运行后输出结果 x = 1.6471e-005这是因为化为后，计算机运行程序 x= (8(-19)/(sqrt(1+8(-19)+1)运行后的结果为x =3.4694e-018由于有效数字的损失甚少，所以运算的结果再与继续计算，最后输出的结果与的精确值相差无几。例1.8 求数的近似值。解 （1）直接用MATLAB程序 x=30;x1= sqrt(x2-1)运行后输出结果x1 = 29.9833输入MATLAB程序 x=30; x1=29.9833;y=log(x-x1)运行后输出结果y = -4.0923（2）因为中的很大，如果采用倒数变换法，即 .输入MATLAB程序 x=30;y=-log(x+sqrt(x2-1)运行后输出结果y = -4.0941（3）输入MATLAB程序 x=30; y=log(x-sqrt(x2-1)运行后输出结果y = -4.0941可见，（2）计算的近似值比（1）的误差小。参加计算的数，有时数量级相差很大.如果不注意采取相应的措施，在它们的加减法运算中，绝对值很小的那个数经常会被绝对值较大的那个数“吃掉”，不能发挥其作用，造成计算结果失真。例1.9 请在16位十进制数值精度计算机上利用软件MATLAB计算下面的两个数和将计算结果与准确值比较，解释计算结果。解 在MATLAB工作窗口输入下面程序 x=111111111111111+0.1+0.3, y=1111111111111111+0.1+0.3运行后输出结果 x = 1.111111111111114e+014，y =1.111111111111111e+015从输出的结果可以看出，x，而y.为什么仅仅比多一位1，而y呢？这是因为计算机进行运算时，首先要把参加运算的数写成绝对值小于1而“阶码”相同的数，这一过程称为数的“对阶”。在16位十进制数值精度计算机上利用软件MATLAB计算这两个数，把运算的数写成浮点规格化形式为在16位十进制数值精度计算机上，三项的数都表示为小数点后面16位数字的数与之积，所以，计算机没有对数进行截断，而是按原来的三个数进行计算。因此，计算的结果。而三项的数都表示写成绝对值小于1而“阶码”都为的数以后，第一项的纯小数的小数点后面有16位数字.但是，后两项数的纯小数的小数点后面有17位数字，超过了16位十进制数值精度计算机的存储量，计算机对后两项的数都进行截断最后一位，即后两项的数都是16位机上的零，再进行计算，所以计算结果与实际不符。三、实验任务对，计算定积分 .算法1：利用递推公式 , ,取 .算法2：利用递推公式 .注意到,取 .思考：从计算结果看，哪个算法是不稳定的，哪个算法是稳定的。实验二、插值法插值法是函数逼近的一种重要方法，它是数值积分、微分方程数值解等数值计算的基础与工具，其中多项式插值是最常用和最基本的方法。拉格朗日插值多项式的优点是表达式简单明确，形式对称，便于记忆，它的缺点是如果想要增加插值节点，公式必须整个改变，这就增加了计算工作量。而牛顿插值多项式对此做了改进，当增加一个节点时只需在原牛顿插值多项式基础上增加一项，此时原有的项无需改变，从而达到节省计算次数、节约存储单元、应用较少节点达到应有精度的目的。一、实验目的1、理解插值的基本概念，掌握各种插值方法，包括拉格朗日插值和牛顿插值等，注意其不同特点；2、通过实验进一步理解并掌握各种插值的基本算法。二、算法实例1. 与插值有关的MATLAB 函数(一) POLY2SYM 函数调用格式一：poly2sym (C)调用格式二：f1=poly2sym(C,V) 或 f2=poly2sym(C, sym (V) ),(二) POLYVAL 函数调用格式：Y = polyval(P,X)(三) POLY 函数调用格式：Y = poly (V)(四) CONV 函数调用格式：C =conv (A, B)(五) DECONV 函数调用格式：Q,R =deconv (B,A)(六) roots(poly(1:n)命令调用格式：roots(poly(1:n) (七) det(a*eye(size (A) - A)命令调用格式：b=det(a*eye(size (A) - A)2.线性插值及其MATLAB程序例2.1 已知函数在上具有二阶连续导数，且满足条件。求线性插值多项式和函数值，并估计其误差。解：输入程序 X=1,3;Y=1,2; l01= poly(X(2)/( X(1)- X(2), l11= poly(X(1)/( X(2)- X(1), l0=poly2sym (l01),l1=poly2sym (l11), P = l01* Y(1)+ l11* Y(2), L=poly2sym (P),x=1.5; Y = polyval(P,x)运行后输出基函数l0和l1及其插值多项式的系数向量P（略）、插值多项式L和插值Y为l0 =-1/2*x+3/2 l1 =1/2*x-1/2 L = 1/2*x+1/2 Y =1.2500输入程序 M=5;R1=M*abs(x-X(1)* (x-X(2)/2运行后输出误差限为 R1 = 1.8750例2.2 求函数e在上线性插值多项式，并估计其误差。解： 输入程序 X=0,1; Y =exp(-X) , l01= poly(X(2)/( X(1)- X(2), l11= poly(X(1)/( X(2)- X(1), l0=poly2sym (l01),l1=poly2sym (l11), P = l01* Y(1)+ l11* Y(2), L=poly2sym (P),运行后输出基函数l0和l1及其插值多项式的系数向量P和插值多项式L为l0 = -x+1 l1 = x P = -0.6321 1.0000L =-1423408956596761/2251799813685248*x+1 输入程序 M=1;x=0:0.001:1; R1=M*max(abs(x-X(1).*(x-X(2)./2运行后输出误差限为 R1 = 0.1250.3.抛物线插值及其MATLAB程序例2.3 求将区间 0, /2 分成等份，用产生个节点，然后分别作线性插值函数和抛物线插值函数.用它们分别计算cos (/6) (取四位有效数字)，并估计其误差。解： 输入程序 X=0,pi/2; Y =cos(X) ,l01= poly(X(2)/( X(1)- X(2), l11= poly(X(1)/( X(2)- X(1), l0=poly2sym (l01),l1=poly2sym (l11), P = l01* Y(1)+ l11* Y(2), L=poly2sym (P),x=pi/6; Y = polyval(P,x)运行后输出基函数l0和l1及其插值多项式的系数向量P、插值多项式和插值为l0 =-5734161139222659/9007199254740992*x+1l1 =5734161139222659/9007199254740992*xP = -0.6366 1.0000L =-5734161139222659/9007199254740992*x+1Y =0.6667输入程序 M=1;x=pi/6; R1=M*abs(x-X(1)*(x-X(2)/2运行后输出误差限为R1 = 0.2742.（2） 输入程序 X=0:pi/4:pi/2; Y =cos(X) ,l01= conv (poly(X(2),poly(X(3)/( X(1)- X(2)* ( X(1)- X(3), l11= conv (poly(X(1), poly(X(3)/( X(2)- X(1)* ( X(2)- X(3),l21= conv (poly(X(1), poly(X(2)/( X(3)- X(1)* ( X(3)- X(2),l0=poly2sym (l01),l1=poly2sym (l11),l2=poly2sym (l21),P = l01* Y(1)+ l11* Y(2) + l21* Y(3), L=poly2sym (P),x=pi/6; Y = polyval(P,x)运行后输出基函数l01、l11和l21及其插值多项式的系数向量P、插值多项式L和插值Y为l0=228155022448185/281474976710656*x2-2150310427208497/1125899906842624*x+1l1=-228155022448185/140737488355328*x2+5734161139222659/2251799813685248*xl2=228155022448185/281474976710656*x2-5734161139222659/9007199254740992*xP = -0.3357 -0.1092 1.0000L=-6048313895780875/18014398509481984*x2-7870612110600739/72057594037927936*x+1Y = 0.8508输入程序 M=1;x=pi/6; R2=M*abs(x-X(1)*(x-X(2) *(x-X(3)/6运行后输出误差限为R2 =0.0239.4.次拉格朗日（Lagrange）插值及其MATLAB程序例2.4 给出节点数据，作三次拉格朗日插值多项式计算，并估计其误差。解：输入程序 X=-2,0,1,2; Y =17,1,2,17;p1=poly(X(1); p2=poly(X(2);p3=poly(X(3); p4=poly(X(4); l01= conv ( conv (p2, p3), p4)/( X(1)- X(2)* ( X(1)- X(3) * ( X(1)- X(4), l11= conv ( conv (p1, p3), p4)/( X(2)- X(1)* ( X(2)- X(3) * ( X(2)- X(4),l21= conv ( conv (p1, p2), p4)/( X(3)- X(1)* ( X(3)- X(2) * ( X(3)- X(4),l31= conv ( conv (p1, p2), p3)/( X(4)- X(1)* ( X(4)- X(2) * ( X(4)- X(3),l0=poly2sym (l01),l1=poly2sym (l11),l2=poly2sym (l21), l3=poly2sym (l31),P = l01* Y(1)+ l11* Y(2) + l21* Y(3) + l31* Y(4),运行后输出基函数l0，l1，l2和l3及其插值多项式的系数向量P（略）为l0 =-1/24*x3+1/8*x2-1/12*x，l1 =1/4*x3-1/4*x2-x+1l2 =-1/3*x3+4/3*x，l3 =1/8*x3+1/8*x2-1/4*x输入程序 L=poly2sym (P),x=0.6; Y = polyval(P,x)运行后输出插值多项式和插值为L = x3+4*x2-4*x+1 Y = 0.2560.输入程序 syms M; x=0.6; R3=M*abs(x-X(1)*(x-X(2) *(x-X(3) *(x-X(4)/24运行后输出误差限为R3 =91/2500*M即 R3 ， .5. 拉格朗日插值及其误差估计的MATLAB程序拉格朗日插值及其误差估计的MATLAB主程序function y,R=lagranzi(X,Y,x,M)n=length(X); m=length(x);for i=1:m z=x(i);s=0.0; for k=1:n p=1.0; q1=1.0; c1=1.0;for j=1:n if j=kp=p*(z-X(j)/(X(k)-X(j); end q1=abs(q1*(z-X(j);c1=c1*j; end s=p*Y(k)+s; end y(i)=s;endR=M*q1/c1;例 2.5 已知，用拉格朗日插值及其误差估计的MATLAB主程序求的近似值，并估计其误差。解：在MATLAB工作窗口输入程序 x=2*pi/9; M=1; X=pi/6 ,pi/4, pi/3;Y=0.5,0.7071,0.8660; y,R=lagranzi(X,Y,x,M)运行后输出插值y及其误差限R为y =0.6434 R =8.8610e-004.6.牛顿插值及其误差估计的MATLAB主程序function y,R= newcz(X,Y,x,M)n=length(X); m=length(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(t)= polyval(C, z);endR=M*q1/c1;例 2.6 已知，用牛顿插值法求的近似值，估计其误差，并与例 6.2.6的计算结果比较.解 方法1（牛顿插值及其误差估计的MATLAB主程序）输入MATLAB程序 x=2*pi/9;M=1; X=pi/6 ,pi/4, pi/3; Y=0.5,0.7071,0.8660; y,R= newcz(X,Y,x,M)运行后输出y = R =0.6434 8.8610e-004方法2 （求牛顿插值多项式和差商的MATLAB主程序）输入MATLAB程序 x=2*pi/9; X=pi/6 ,pi/4, pi/3; Y=0.5,0.7071,0.8660; M=1; A,C,L,wcgs,Cw= newploy(X,Y), y=polyval(C,x)运行后输出结果A = 0.5000 0 0 0.7071 0.7911 0 0.8660 0.6070 -0.3516C = -0.3516 1.2513 -0.0588L =-1583578379808357/4503599627370496*x2+1408883391907005/1125899906842624*x-132405829044691/2251799813685248wcgs =1/6*M*(x3-3/4*x2*pi+4012734077357799/2251799813685248*x-7757769783530263/18014398509481984)Cw = 0.1667 -0.3927 0.2970 -0.0718y = 0.6434上述两种方法计算y的结果相同.7. 牛顿插值法的MATLAB综合程序求牛顿插值多项式、差商、插值及其误差估计的MATLAB主程序function y,R,A,C,L=newdscg(X,Y,x,M)n=length(X); m=length(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(t)= polyval(C, z);endR=M*q1/c1;L(k,:)=poly2sym(C);例 2.7 给出节点数据，作三阶牛顿插值多项式，计算，并估计其误差。解 首先将名为newdscg.m的程序保存为M文件，然后在MATLAB工作窗口输入程序 syms M,X=-4,0,1,2; Y =27,1,2,17; x=-2.345; y,R,A,C,P=newdscg(X,Y,x,M)运行后输出插值y及其误差限公式R，三阶牛顿插值多项式P及其系数向量C，差商的矩阵A如下y = 22.3211R =1323077530165133/562949953421312*M（即R =2.3503*M）A= 27.0000 0 0 0 1.0000 -6.5000 0 0 2.0000 1.0000 1.5000 0 17.0000 15.0000 7.0000 0.9167C =0.9167 4.2500 -4.1667 1.0000P =11/12*x3+17/4*x2-25/6*x+1例 2.8 将区间 0,/2 分成等份，用产生个节点，求二阶和三阶牛顿插值多项式和.用它们分别计算sin (/7) (取四位有效数字)，并估计其误差.解 首先将名为newdscg.m的程序保存为M文件，然后在MATLAB工作窗口输入程序 X1=0:pi/4:pi/2; Y1 =sin(X1); M=1; x=pi/7; X2=0:pi/6:pi/2; Y2 =sin(X2); y1,R1,A1,C1,P2=newdscg(X1,Y1,x,M), y2,R2,A2,C2,P3=newdscg(X2,Y2,x,M)运行后输出插值y1=和y2=及其误差限R1和R2，二阶和三阶牛顿插值多项式P2和P3及其系数向量C1和C2，差商的矩阵A1和A2如下y1 =0.4548R1 =0.0282A1 = 0 0 0 0.7071 0.9003 0 1.0000 0.3729 -0.3357C1 = -0.3357 1.1640 0P2=-3024156947890437/9007199254740992*x2+163820246322191/140737488355328*xy2 =0.4345R2 =9.3913e-004A2 = 0 0 0 0 0.5000 0.9549 0 0 0.8660 0.6991 -0.2443 0 1.0000 0.2559 -0.4232 -0.1139C2 =-0.1139 -0.0655 1.0204 0P3=-1025666884451963/9007199254740992*x3-4717668559161127/72057594037927936*x2+4595602396951547/4503599627370496*x三、实验任务1、 已知函数表 0.56160 0.56280 0.56401 0.56521 0.82741 0.82659 0.82577 0.82495用二次拉格朗日插值多项式求时的函数近似值。2、 已知函数表 0.4 0.55 0.65 0.8 0.9 0.41075 0.57815 0.69675 0.88811 1.02652用牛顿插值多项式求和。四、思考题1、 插值多项式的存在唯一性有何意义？2、 多项式插值是如何构造的？截断误差如何表示？如何估计？3、 在插值区间上，随着节点的增多，插值多项式是否越来越接近被插函数？实验三、数据拟合法曲线拟合的最小二乘法是计算机数据处理的重要内容，也是函数逼近的另一种重要方法，它在工程技术中有着广泛的应用。对实际问题而言，拟合曲线的选择是一个极其重要而又比较困难的问题，必要时可由草图观察选取几种不同类型的拟合曲线，再以其偏差小者为优，经检验后再决定最后的取舍。一、实验目的1、 理解数据拟合的基本概念、基本方法；2、 掌握最小二乘法的基本原理，并会通过计算机解决实际问题；3、了解超定方程组的最小二乘解法。二、算法实例例3.1 给出一组数据点列入表31中，试用线性最小二乘法求拟合曲线，并估计其误差，作出拟合曲线。表31 例3.1的一组数据xi-2.5 -1.7 -1.1 -0.8 0 0.1 1.5 2.7 3.6yi-192.9 -85.50 -36.15 -26.52 -9.10 -8.43 -13.12 6.50 68.04解 （1）在MATLAB工作窗口输入程序 x=-2.5 -1.7 -1.1 -0.8 0 0.1 1.5 2.7 3.6; y=-192.9 -85.50 -36.15 -26.52 -9.10 -8.43 -13.12 6.50 68.04;plot(x,y,r*),legend(实验数据(xi,yi)xlabel(x), ylabel(y),title(例3.1的数据点(xi,yi)的散点图)运行后屏幕显示数据的散点图（略）。 （3）编写下列MATLAB程序计算在处的函数值，即输入程序 syms a1 a2 a3 a4x=-2.5 -1.7 -1.1 -0.8 0 0.1 1.5 2.7 3.6; fi=a1.*x.3+ a2.*x.2+ a3.*x+ a4运行后屏幕显示关于a1,a2, a3和a4的线性方程组fi = -125/8*a1+25/4*a2-5/2*a3+a4, -4913/1000*a1+289/100*a2-17/10*a3+a4, -1331/1000*a1+121/100*a2-11/10*a3+a4, -64/125*a1+16/25*a2-4/5*a3+a4, a4,1/1000*a1+1/100*a2+1/10*a3+a4, 27/8*a1+9/4*a2+3/2*a3+a4, 19683/1000*a1+729/100*a2+27/10*a3+a4, 5832/125*a1+324/25*a2+18/5*a3+a4编写构造误差平方和的MATLAB程序 y=-192.9 -85.50 -36.15 -26.52 -9.10 -8.43 -13.12 6.50 68.04;fi=-125/8*a1+25/4*a2-5/2*a3+a4, -4913/1000*a1+289/100*a2-17/10*a3+a4, -1331/1000*a1+121/100*a2-11/10*a3+a4, -64/125*a1+16/25*a2-4/5*a3+a4,a4, 1/1000*a1+1/100*a2+1/10*a3+a4, 27/8*a1+9/4*a2+3/2*a3+a4, 19683/1000*a1+729/100*a2+27/10*a3+a4, 5832/125*a1+324/25*a2+18/5*a3+a4;fy=fi-y; fy2=fy.2; J=sum(fy.2)运行后屏幕显示误差平方和如下J=(-125/8*a1+25/4*a2-5/2*a3+a4+1929/10)2+(-4913/1000*a1+289/100*a2-17/10*a3+a4+171/2)2+(-1331/1000*a1+121/100*a2-11/10*a3+a4+723/20)2+(-64/125*a1+16/25*a2-4/5*a3+a4+663/25)2+(a4+91/10)2+(1/1000*a1+1/100*a2+1/10*a3+a4+843/100)2+(27/8*a1+9/4*a2+3/2*a3+a4+328/25)2+(19683/1000*a1+729/100*a2+27/10*a3+a4-13/2)2+(5832/125*a1+324/25*a2+18/5*a3+a4-1701/25)2为求使达到最小，只需利用极值的必要条件 ，得到关于的线性方程组，这可以由下面的MATLAB程序完成，即输入程序 syms a1 a2 a3 a4J=(-125/8*a1+25/4*a2-5/2*a3+a4+1929/10)2+(-4913/1000*a1+289/100*a2-17/10*a3+a4.+171/2)2+(-1331/1000*a1+121/100*a2-11/10*a3+a4+723/20)2+(-64/125*a1+16/25*a2-4/5*a3+a4+663/25)2+(a4+91/10)2+(1/1000*a1+1/100*a2+1/10*a3+a4+843/100)2+(27/8*a1+9/4*a2+3/2*a3+a4+328/25)2+(19683/1000*a1+729/100*a2+27/10*a3+a4-13/2)2+(5832/125*a1+324/25*a2+18/5*a3+a4-1701/25)2; Ja1=diff(J,a1); Ja2=diff(J,a2); Ja3=diff(J,a3); Ja4=diff(J,a4);Ja11=simple(Ja1), Ja21=simple(Ja2), Ja31=simple(Ja3), Ja41=simple(Ja4),运行后屏幕显示J分别对a1, a2 ,a3 ,a4的偏导数如下Ja11=56918107/10000*a1+32097579/25000*a2+1377283/2500*a3+23667/250*a4-8442429/625Ja21=32097579/25000*a1+1377283/2500*a2+23667/250*a3+67*a4+767319/625Ja31=1377283/2500*a1+23667/250*a2+67*a3+18/5*a4-232638/125Ja41=23667/250*a1+67*a2+18/5*a3+18*a4+14859/25解线性方程组Ja11 =0，Ja21 =0，