数值分析上机实验报告

实验 报告 软件系统测试报告下载sgs报告如何下载关于路面塌陷情况报告535n,sgs报告怎么下载竣工报告下载 yyhhit@163.com .1. 实验报告一 题目： 非线性方程求解 摘要：非线性方程的解析解通常很难给出，因此线性方程的数值解法就尤为重要。本实验 采用两种常见的求解方法二分法和 Newton 法及改进的 Newton 法。 前言：（目的和意义） 掌握二分法与 Newton 法的基本原理和应用。 数学原理： 对于一个非线性方程的数值解法很多。在此介绍两种最常见的方法：二分法和 Newton 法。 对于二分法，其数学实质就是说对于给定的待求解的方程f(x)，其在[a,b]上连续， f(a)f(b)<0，且f(x)在[a,b]内仅有一个实根x*，取区间中点c，若，则c恰为其根，否则根据 f(a)f(c)<0 是否成立判断根在区间[a,c]和[c,b]中的哪一个，从而得出新区间，仍称为[a,b]。 重复运行计算，直至满足精度为止。这就是二分法的计算思想。 Newton法通常预先要给出一个猜测初值x0，然后根据其迭代公式 )( )( '1 k k kk xf xfxx −=+ 产生逼近解x*的迭代数列{xk}，这就是Newton法的思想。当x0接近x*时收敛很快，但是当 x0选择不好时，可能会发散，因此初值的选取很重要。另外，若将该迭代公式改进为 )( )( '1 k k kk xf xfrxx −=+ 其中 r 为要求的方程的根的重数，这就是改进的 Newton 法，当求解已知重数的方程的根 时，在同种条件下其收敛速度要比 Newton 法快的多。 程序 设计 领导形象设计圆作业设计ao工艺污水处理厂设计附属工程施工组织设计清扫机器人结构设计 ： 本实验采用 Matlab 的 M 文件编写。其中待求解的方程写成 function 的方式，如下 function y=f(x); y=-x*x-sin(x); 写成如上形式即可，下面给出主程序。 二分法源程序： clear 实验报告 yyhhit@163.com .2. %%%给定求解区间 b=1.5; a=0; %%%误差 R=1; k=0;%迭代次数初值 while (R>5e-6) ; c=(a+b)/2; if f12(a)*f12(c)>0; a=c; else b=c; end R=b-a;%求出误差 k=k+1; end x=c%给出解 Newton 法及改进的 Newton 法源程序： clear %%%% 输入函数 f=input('请输入需要求解函数>>','s') %%%求解 f(x)的导数 df=diff(f); %%%改进常数或重根数 miu=2; %%%初始值 x0 x0=input('input initial value x0>>'); k=0;%迭代次数 max=100;%最大迭代次数 R=eval(subs(f,'x0','x'));%求解 f(x0)，以确定初值 x0 时否就是解 while (abs(R)>1e-8) x1=x0-miu*eval(subs(f,'x0','x'))/eval(subs(df,'x0','x')); R=x1-x0; x0=x1; k=k+1; if (eval(subs(f,'x0','x'))<1e-10); 实验报告 yyhhit@163.com .3. break end if k>max;%如果迭代次数大于给定值，认为迭代不收敛，重新输入初值 ss=input('maybe result is error,choose a new x0,y/n?>>','s'); if strcmp(ss,'y') x0=input('input initial value x0>>'); k=0; else break end end end k;%给出迭代次数 x=x0；%给出解 结果分析和讨论： 1. 用二分法计算方程 0 2 sin 2 =− xx 在[1，2]内的根。( ,下同) 610*5 −=ε 计算结果为 x= 1.40441513061523； f(x)= -3.797205105904311e-007； k=18； 由 f(x)知结果满足要求，但迭代次数比较多，方法收敛速度比较慢。 2. 用二分法计算方程 在[1，1.5]内的根。 013 =−− xx 计算结果为 x= 1.32471847534180； f(x)= 2.209494846194815e-006； k=17； 由 f(x)知结果满足要求，但迭代次数还是比较多。 3. 用 Newton 法求解下列方程 a) x01 =−xxe 0=0.5； 计算结果为 x= 0.56714329040978； f(x)= 2.220446049250313e-016； k=4； 由 f(x)知结果满足要求，而且又迭代次数只有 4 次看出收敛速度很快。 实验报告 yyhhit@163.com .4. b) x013 =−− xx 0=1； c) x0)12()1( 2 =−− xx 0=0.45, x0=0.65； 当x0=0.45 时，计算结果为 x= 0.49999999999983； f(x)= -8.362754932994584e-014； k=4； 由 f(x)知结果满足要求，而且又迭代次数只有 4 次看出收敛速度很快，实际上该方程确实 有真解 x=0.5。 当x0=0.65 时，计算结果为 x= 0.50000000000000； f(x)=0； k=9； 由f(x)知结果满足要求，实际上该方程确实有真解x=0.5，但迭代次数增多，实际上当取x0〉 0.68 时，x≈1，就变成了方程的另一个解，这说明Newton法收敛与初值很有关系，有的时 候甚至可能不收敛。 4. 用改进的 Newton 法求解，有 2 重根，取 2=µ 0)12()1( 2 =−− xx x0=0.55；并与 3.中的c)比较结果。 当x0=0.55 时，程序死循环，无法计算，也就是说不收敛。改 5.1=µ 时，结果收敛为 x=0.50000087704286； f(x)=4.385198907621127e-007； k=16； 显然这个结果不是很好，而且也不是收敛至方程的 2 重根上。 当x0=0.85 时，结果收敛为 x= 1.00000000000489； f(x)= 2.394337647718737e-023； k=4； 这次达到了预期的结果，这说明初值的选取很重要，直接关系到方法的收敛性，实际上直 接用 Newton 法，在给定同样的条件和精度要求下，可得其迭代次数 k=15，这说明改进后 的 Newton 法法速度确实比较快。 结论： 对于二分法，只要能够保证在给定的区间内有根，使能够收敛的，当时收敛的速度和 给定的区间有关，二且总体上来说速度比较慢。Newton 法，收敛速度要比二分法快，但 是最终其收敛的结果与初值的选取有关，初值不同，收敛的结果也可能不一样，也就是结 果可能不时预期需要得结果。改进的 Newton 法求解重根问题时，如果初值不当，可能会 不收敛，这一点非常重要，当然初值合适，相同情况下其速度要比 Newton 法快得多。 实验报告 yyhhit@163.com .5. 实验报告二 题目： Gauss 列主元消去法 摘要：求解线性方程组的方法很多，主要分为直接法和间接法。本实验运用直接法的 Guass 消去法,并采用选主元的方法对方程组进行求解。 前言：（目的和意义） 1. 学习 Gauss 消去法的原理。 2. 了解列主元的意义。 3. 确定什么时候系数阵要选主元 数学原理： 由于一般线性方程在使用 Gauss 消去法求解时，从求解的过程中可以看到，若 =0， 则必须进行行交换，才能使消去过程进行下去。有的时候即使 0，但是其绝对值非 常小，由于机器舍入误差的影响，消去过程也会出现不稳定得现象，导致结果不正确。因 此有必要进行列主元技术，以最大可能的消除这种现象。这一技术要寻找行 r，使得 )1( −k kka ≠− )1(kkka )1()1( max|| − > − = kikki k rk aa 并将第 r 行和第 k 行的元素进行交换，以使得当前的 的数值比 0 要大的多。这种列主 元的消去法的主要步骤如下： )1( −k kka 1. 消元过程 对 k=1,2,…,n-1,进行如下步骤。 1) 选主元，记 ikkirk aa > = max|| 若 很小，这说明方程的系数矩阵严重病态，给出警告，提示结果可能不对。 || rka 2) 交换增广阵 A 的 r，k 两行的元素。 kjrj aa ↔ (j=k,…,n+1) 3) 计算消元 kkkjikijij aaaaa /−= (i=k+1,…,n; j=k+1,……,n+1) 2. 回代过程 对 k= n, n-1,…,1,进行如下计算 实验报告 yyhhit@163.com .6. )/( 1 1, ∑ −= + −= n kj kkjkjnkk axaax 至此，完成了整个方程组的求解。 程序设计： 本实验采用 Matlab 的 M 文件编写。 Gauss 消去法源程序： clear a=input('输入系数阵：>>\n') b=input('输入列阵 b：>>\n') n=length(b); A=[a b] x=zeros(n,1); %%%函数主体 for k=1:n-1; %%%是否进行主元选取 if abs(A(k,k))abs(t) p=r; else p=k; end end %%%交换元素 if p~=k; for q=k:n+1; s=A(k,q); A(k,q)=A(p,q); A(p,q)=s; end 实验报告 yyhhit@163.com .7. end end %%%判断系数矩阵是否奇异或病态非常严重 if abs(A(k,k))< yipusilong disp(‘矩阵奇异，解可能不正确’) end %%%%计算消元,得三角阵 for r=k+1:n; m=A(r,k)/A(k,k); for q=k:n+1; A(r,q)=A(r,q)-A(k,q)*m; end end end %%%%求解 x x(n)=A(n,n+1)/A(n,n); for k=n-1:-1:1; s=0; for r=k+1:n; s=s+A(k,r)*x(r); end t=(A(k,n+1)-s) x(k)=(A(k,n+1)-s)/A(k,k) end 结果分析和讨论： 例：求解方程 。其中 ⎥⎥ ⎥ ⎦ ⎤ ⎢⎢ ⎢ ⎣ ⎡ = ⎥⎥ ⎥ ⎦ ⎤ ⎢⎢ ⎢ ⎣ ⎡ ⎥⎥ ⎥ ⎦ ⎤ ⎢⎢ ⎢ ⎣ ⎡ 10 34 22 123 575 62 z y xε ε 为一小数，当 时， 分别采用列主元和不列主元的 Gauss 消去法求解，并比较结果。 2014105 10,10,10,10 −−−−=ε 记Emax为求出的解代入方程后的最大误差，按要求，计算结果如下： 当 时，不选主元和选主元的计算结果如下，其中前一列为不选主元结果，后 一列为选主元结果，下同。 510−=ε 0.99999934768391 0.99999934782651 2.00000217421972 2.00000217391163 2.99999760859451 2.99999760869721 实验报告 yyhhit@163.com .8. Emax= 9.301857062382624e-010，0 此时，由于ε 不是很小，机器误差就不是很大，由Emax可以看出不选主元的计算结果 精度还可以，因此此时可以考虑不选主元以减少计算量。 当 时，不选主元和选主元的计算结果如下 1010−=ε 1.00001784630877 0.99999999999348 1.99998009720807 2.00000000002174 3.00000663424731 2.99999999997609 Emax= 2.036758973744668e-005，0 此时由Emax可以看出不选主元的计算精度就不好了，误差开始增大。 当 时，不选主元和选主元的计算结果如下 1410−=ε 1.42108547152020 1.00000000000000 1.66666666666666 2.00000000000000 3.11111111111111 300000000000000 Emax= 0.70770085900503，0 此时由Emax可以看出，不选主元的结果应该可以说是不正确了，这是由机器误差引起 的。 当 时，不选主元和选主元的计算结果如下 2010−=ε NaN 1 NaN 2 NaN 3 Emax=NaN, 0 不选主元时，程序报错：Warning: Divide by zero.。这是因为机器计算的最小精度为 10-15，所以此时的 就认为是 0，故出现了错误现象。而选主元时则没有这种现象， 而且由E 2010−=ε max可以看出选主元时的结果应该是精确解。 结论： 采用Gauss消去法时，如果在消元时对角线上的元素始终较大（假如大于 10-5），那么 本方法不需要进行列主元计算，计算结果一般就可以达到要求，否则必须进行列主元这一 步，以减少机器误差带来的影响，使方法得出的结果正确。 实验报告 yyhhit@163.com .9. 实验报告三 题目： Rung 现象产生和克服 摘要：由于高次多项式插值不收敛，会产生 Runge 现象，本实验在给出具体的实例后， 采用分段线性插值和三次样条插值的方法有效的克服了这一现象，而且还取的很好的插值 效果。 前言：（目的和意义） 1. 深刻认识多项式插值的缺点。 2. 明确插值的不收敛性怎样克服。 3. 明确精度与节点和插值方法的关系。 数学原理： 在给定 n+1 个节点和相应的函数值以后构造 n 次的 Lagrange 插值多项式，实验结果 表明（见后面的图）这种多项式并不是随着次数的升高对函数的逼近越来越好，这种现象 就是 Rung 现象。 解决 Rung 现象的方法通常有分段线性插值、三次样条插值等方法。 分段线性插值： 设在区间[a, b]上，给定 n+1 个插值节点 a=x0s(n); l=(x-s(n))./(s(n+1)-s(n)); else l=0; end otherwise if x>=s(k-1)&x<=s(k); l=(x-s(k-1))./(s(k)-s(k-1)); else if x>=s(k)&x<=s(k+1); l=(x-s(k+1))./(s(k)-s(k+1)); else l=0; end end end ff=ff+Rf(s(k))*l;%%求插值函数值 end m=m+1; f(m)=ff; end %%%作出曲线 x=-1:hh:1; plot(x,f,'r'); grid on hold on 三次样条插值源程序：（采用第一边界条件） clear n=input('将区间分为的等份数输入：\n'); %%%插值区间 a=-1; b=1; hh=0.001;%画图的步长 s=[a+(b-a)/n*[0:n]];%%%给定的定点，Rf 为给定的函数 %%%%第一边界条件 Rf"(-1),Rf"(1) 实验报告 yyhhit@163.com .12. v=5000*1/(1+25*a*a)^3-50/(1+25*a*a)^4; for k=1:n;%取出节点间距 h(k)=s(k+1)-s(k); end for k=1:n-1;%求出系数向量 lamuda,miu la(k)=h(k+1)/(h(k+1)+h(k)); miu(k)=1-la(k); end %%%%赋值系数矩阵 A for k=1:n-1; for p=1:n-1; switch p case k A(k,p)=2; case k-1 A(k,p)=miu(p+1); case k+1 A(k,p)=la(p-1); otherwise A(k,p)=0; end end end %%%%求出 d 阵 for k=1:n-1; switch k case 1 d(k)=6*f2c([s(k) s(k+1) s(k+2)])-miu(k)*v; case n-1 d(k)=6*f2c([s(k) s(k+1) s(k+2)])-la(k)*v; otherwise d(k)=6*f2c([s(k) s(k+1) s(k+2)]); end end %%%%求解 M 阵 M=A\d'; 实验报告 yyhhit@163.com .13. M=[v;M;v]; %%%% m=0; f=0; for x=a:hh:b; if x==a; p=1; else p=ceil((x-s(1))/((b-a)/n)); end ff1=0; ff2=0; ff3=0; ff4=0; m=m+1; ff1=1/h(p)*(s(p+1)-x)^3*M(p)/6; ff2=1/h(p)*(x-s(p))^3*M(p+1)/6; ff3=((Rf(s(p+1))-Rf(s(p)))/h(p)-h(p)*(M(p+1)-M(p))/6)*(x-s(p)); ff4=Rf(s(p))-M(p)*h(p)*h(p)/6; f(m)=ff1+ff2+ff3+ff4 ; end %%%作出插值图形 x=a:hh:b; plot(x,f,'k') hold on grid on 结果分析和讨论： 本实验采用函数 2251 1)( x xf += 进行数值插值，插值区间为[-1，1]，给定节点为 xj=-1+jh,h=0.1，j=0,…，n。下面分别给出Lagrange插值，三次样条插值，线性插值的函 数曲线和数据表。图中只标出Lagrange插值的十次多项式的曲线，其它曲线没有标出，从 数据表中可以看出具体的误差。 实验报告 yyhhit@163.com .14. 表中，L10(x)为Lagrange插值的 10 次多项式，S10(x)，S40(x)分别代表n=10，40 的三次样条 插值函数，X10(x)，X40(x)分别代表n=10，40 的线性分段插值函数。 x f(x) L10(x) S10(x) S40(x) X10(x) X40(x) -1.00000000000000 0.03846153846154 0.03846153846154 0.03846153846154 0.03846153846154 0.03846153846154 0.03846153846154 -0.95000000000000 0.04244031830239 1.92363114971920 0.04240833151040 0.04244031830239 0.04355203619910 0.04244031830239 -0.90000000000000 0.04705882352941 1.57872099034926 0.04709697585458 0.04705882352941 0.04864253393665 0.04705882352941 -0.85000000000000 0.05245901639344 0.71945912837982 0.05255839923979 0.05245901639344 0.05373303167421 0.05245901639344 -0.80000000000000 0.05882352941176 0.05882352941176 0.05882352941176 0.05882352941176 0.05882352941176 0.05882352941176 -0.75000000000000 0.06639004149378 -0.23146174989674 0.06603986172744 0.06639004149378 0.06911764705882 0.06639004149378 -0.70000000000000 0.07547169811321 -0.22619628906250 0.07482116198866 0.07547169811321 0.07941176470588 0.07547169811321 -0.65000000000000 0.08648648648649 -0.07260420322418 0.08589776360849 0.08648648648649 0.08970588235294 0.08648648648649 -0.60000000000000 0.10000000000000 0.10000000000000 0.10000000000000 0.10000000000000 0.10000000000000 0.10000000000000 -0.55000000000000 0.11678832116788 0.21559187891257 0.11783833017713 0.11678832116788 0.12500000000000 0.11678832116788 -0.50000000000000 0.13793103448276 0.25375545726103 0.14004371555730 0.13793103448276 0.15000000000000 0.13793103448276 -0.45000000000000 0.16494845360825 0.23496854305267 0.16722724315883 0.16494845360825 0.17500000000000 0.16494845360825 -0.40000000000000 0.20000000000000 0.20000000000000 0.20000000000000 0.20000000000000 0.20000000000000 0.20000000000000 -0.35000000000000 0.24615384615385 0.19058046675376 0.24054799403464 0.24615384615385 0.27500000000000 0.24615384615385 -0.30000000000000 0.30769230769231 0.23534659131080 0.29735691695860 0.30769230769231 0.35000000000000 0.30769230769231 -0.25000000000000 0.39024390243902 0.34264123439789 0.38048738140327 0.39024390243902 0.42500000000000 0.39024390243902 -0.20000000000000 0.50000000000000 0.50000000000000 0.50000000000000 0.50000000000000 0.50000000000000 0.50000000000000 -0.15000000000000 0.64000000000000 0.67898957729340 0.65746969368431 0.64000000000000 0.62500000000000 0.64000000000000 -0.10000000000000 0.80000000000000 0.84340742982890 0.82052861660828 0.80000000000000 0.75000000000000 0.80000000000000 -0.05000000000000 0.94117647058824 0.95862704866073 0.94832323122810 0.94117647058824 0.87500000000000 0.94117647058824 0 1.00000000000000 1.00000000000000 1.00000000000000 1.00000000000000 1.00000000000000 1.00000000000000 0.05000000000000 0.94117647058824 0.95862704866073 0.94832323122810 0.94117647058824 0.87500000000000 0.94117647058824 0.10000000000000 0.80000000000000 0.84340742982890 0.82052861660828 0.80000000000000 0.75000000000000 0.80000000000000 0.15000000000000 0.64000000000000 0.67898957729340 0.65746969368431 0.64000000000000 0.62500000000000 0.64000000000000 0.20000000000000 0.50000000000000 0.50000000000000 0.50000000000000 0.50000000000000 0.50000000000000 0.50000000000000 0.25000000000000 0.39024390243902 0.34264123439789 0.38048738140327 0.39024390243902 0.42500000000000 0.39024390243902 0.30000000000000 0.30769230769231 0.23534659131080 0.29735691695860 0.30769230769231 0.35000000000000 0.30769230769231 0.35000000000000 0.24615384615385 0.19058046675376 0.24054799403464 0.24615384615385 0.27500000000000 0.24615384615385 0.40000000000000 0.20000000000000 0.20000000000000 0.20000000000000 0.20000000000000 0.20000000000000 0.20000000000000 0.45000000000000 0.16494845360825 0.23496854305267 0.16722724315883 0.16494845360825 0.17500000000000 0.16494845360825 0.50000000000000 0.13793103448276 0.25375545726103 0.14004371555730 0.13793103448276 0.15000000000000 0.13793103448276 实验报告 yyhhit@163.com .15. 0.55000000000000 0.11678832116788 0.21559187891257 0.11783833017713 0.11678832116788 0.12500000000000 0.11678832116788 0.60000000000000 0.10000000000000 0.10000000000000 0.10000000000000 0.10000000000000 0.10000000000000 0.10000000000000 0.65000000000000 0.08648648648649 -0.07260420322418 0.08589776360849 0.08648648648649 0.08970588235294 0.08648648648649 0.70000000000000 0.07547169811321 -0.22619628906250 0.07482116198866 0.07547169811321 0.07941176470588 0.07547169811321 0.75000000000000 0.06639004149378 -0.23146174989674 0.06603986172744 0.06639004149378 0.06911764705882 0.06639004149378 0.80000000000000 0.05882352941176 0.05882352941176 0.05882352941176 0.05882352941176 0.05882352941176 0.05882352941176 0.85000000000000 0.05245901639344 0.71945912837982 0.05255839923979 0.05245901639344 0.05373303167421 0.05245901639344 0.90000000000000 0.04705882352941 1.57872099034926 0.04709697585458 0.04705882352941 0.04864253393665 0.04705882352941 0.95000000000000 0.04244031830239 1.92363114971920 0.04240833151040 0.04244031830239 0.04355203619910 0.04244031830239 1.00000000000000 0.03846153846154 0.03846153846154 0.03846153846154 0.03846153846154 0.03846153846154 0.03846153846154 从以上结果可以看到，用三次样条插值和线性分段插值，不会出现多项式插值是出现 的 Runge 现象，插值效果明显提高。进一步说，为了提高插值精度，用三次样条插值和 线性分段插值是可以增加插值节点的办法来满足要求，而用多项式插值函数时，节点数的 增加必然会使多项式的次数增加，这样会引起数值不稳定，所以说这两种插值要比多项式 插值好的多。而且在给定节点数的条件下，三次样条插值的精度要优于线性分段插值，曲 线的光滑性也要好一些。 实验报告 yyhhit@163.com .16. 实验报告四 题目： 多项式最小二乘法 摘要：对于具体实验时，通常不是先给出函数的解析式，再进行实验，而是通过实验的观 察和测量给出离散的一些点，再来求出具体的函数解析式。又因为测量误差的存在，实际 真实的解析式曲线并不一定通过测量给出的所有点。最小二乘法是求解这一问题的很好的 方法，本实验运用这一方法实现对给定数据的拟合。 前言：（目的和意义） 1. 学习使用最小二成法的原理 2. 了解法方程的特性 数学原理： 对于给定的测量数据(xi,fi)(i=1,2,…，n)，设函数分布为 ∑ = = m j jj xaxy 0 )()( ϕ 特别的，取 )(xjϕ 为多项式 j j xx =)(ϕ (j=0, 1,…，m) 则根据最小二乘法原理，可以构造泛函 ∑ ∑ = = −= n i m j ijjim xafaaaH 1 0 10 ))((),,,( ϕL 令 0=∂ ∂ ka H (k=0, 1,…，m) 则可以得到法方程 ⎥⎥ ⎥⎥ ⎦ ⎤ ⎢⎢ ⎢⎢ ⎣ ⎡ = ⎥⎥ ⎥⎥ ⎦ ⎤ ⎢⎢ ⎢⎢ ⎣ ⎡ ⎥⎥ ⎥⎥ ⎦ ⎤ ⎢⎢ ⎢⎢ ⎣ ⎡ ),( ),( ),( ),(),(),( ),(),(),( ),(),(),( 1 0 1 0 10 11110 00100 mmmmmm m m f f f a a a ϕ ϕ ϕ ϕϕϕϕϕϕ ϕϕϕϕϕϕ ϕϕϕϕϕϕ MM L MMM L L 求该解方程组，则可以得到解 ，因此可得到数据的最小二乘解 maaa ,,, 10 L ∑ = ≈ m j jj xaxf 0 )()( ϕ 实验报告 yyhhit@163.com .17. 程序设计： 本实验采用 Matlab 的 M 文件编写。其中多项式函数 写成 function 的方式，如 下 j j x=ϕ function y=fai(x,j) y=1; for i=1:j y=x.*y; end 写成如上形式即可，下面给出主程序。 多项式最小二乘法源程序 clear %%%给定测量数据点(s,f) s=[3 4 5 6 7 8 9]; f=[2.01 2.98 3.50 5.02 5.47 6.02 7.05]; %%%计算给定的数据点的数目 n=length(f); %%%给定需要拟合的数据的最高次多项式的次数 m=10; %%%程序主体 for k=0:m; g=zeros(1,m+1); for j=0:m; t=0; for i=1:n;%计算内积(fai(si),fai(si)) t=t+fai(s(i),j)*fai(s(i),k); end g(j+1)=t; end A(k+1,:)=g;%法方程的系数矩阵 t=0; for i=1:n;%计算内积(f(si),fai(si)) t=t+f(i)*fai(s(i),k); end b(k+1,1)=t; end 实验报告 yyhhit@163.com .18. a=A\b%求出多项式系数 x=[s(1):0.01:s(n)]'; y=0; for i=0:m; y=y+a(i+1)*fai(x,i); end plot(x,y)%作出拟合成的多项式的曲线 grid on hold on plot(s,f,'rx') %在上图中标记给定的点 结果分析和讨论： 例 用最小二乘法处理下面的实验数据. xi 3 4 5 6 7 8 9 fi 2.01 2.98 3.50 5.02 5.47 6.02 7.05 并作出 的近似分布图。 )(xf 分别采用一次，二次和五次多项式来拟合数据得到相应的拟合多项式为： y1=-0.38643+0.82750x； y2=-1.03024+1.06893x-0.02012x2； y5=-50.75309+51.53527x-19.65947x2+3.66585x3-0.32886x4+0.01137x5； 分别作出它们的曲线图，图中点划线为 y1 曲线，实线为 y2 曲线，虚线为 y5 曲线。’x’为 给定的数据点。从图中可以看出并不是多项式次数越高越好，次数高了，曲线越能给定点 处和实际吻合，但别的地方就很差了。因此，本例选用一次和两次的多项式拟合应该就可 以了。 实验报告 yyhhit@163.com .19. 实验报告五 题目： Romberg 积分法 摘要：对于实际的工程积分问题，很难应用 Newton-Leibnitz 公式去求解。因此应用数值 方法进行求解积分问题已经有着很广泛的应用，本文基于 Romberg 积分法来解决一类积 分问题。 前言：（目的和意义） 1. 理解和掌握 Romberg 积分法的原理； 2. 学会使用 Romberg 积分法； 3. 明确 Romberg 积分法的收敛速度及应用时容易出现的问题。 数学原理： 考虑积分 ，欲求其近似值，通常有复化的梯形公式、Simpsion公式和 Cotes公式。但是给定一个精度，这些公式达到要求的速度很缓慢。如何提高收敛速度， 自然是人们极为关心的课题。为此，记T ∫= ba dxxffI )()( 1,k为将区间[a,b]进行 2k等分的复化的梯形公式计 算结果，记T2,k为将区间[a,b]进行 2k等分的复化的Simpsion公式计算结果，记T3,k为将区间 [a,b]进行 2k等分的复化的Cotes公式计算结果。根据Richardson外推加速方法，可以得到收 敛速度较快的Romberg积分法。其具体的计算公式为： 1. 准备初值，计算 )]()([ 21,1 bfafbaT +−= 2. 按梯形公式的递推关系，计算 ∑− = −+ − +−+−+= 12 0 1,11,1 1 ))5.0( 2 ( 22 1 k i kkkk i abafabTT 3. 按 Romberg 积分公式计算加速值 14 4 1 ,11,1 1 , − −= − −−−+− − − m mkmmkm m mkm TT T m=2,…，k 4. 精度控制。对给定的精度 R ，若 RTT mm <− − 1,11, 则终止计算，并取 为所求结果；否则返回 2 重复计算，直至满足要求的精度为止。 1,mT 实验报告 yyhhit@163.com .20. 程序设计： 本实验采用 Matlab 的 M 文件编写。其中待积分的函数写成 function 的方式，例如如 下 function yy=f(x,y); yy=x.^3; 写成如上形式即可，下面给出主程序 Romberg 积分法源程