Romberg数值积分实验报告班级：__ 姓名：学号：日期：_一、实验目的1、加深外推法的原理理解, 掌握Romberg外推法的计算方法。
2、用matlab软件实现Romberg数值积分来计算题目的运算。
二、基本理论及背景1、理论推导：对区间[a，b]，令h=b-a构造梯形值序列{T2K}。
T1=h[f(a)+f(b)]/2把区间二等分，每个小区间长度为h/2=(b-a)/2，于是T2 =T1/2+[h/2]f(a+h/2）把区间四(2)等分，每个小区间长度为h/2 =（b-a）/4，于是T4 =T2/2+[h/2][f（a+h/4）+f（a+3h/4).....................把[a，b] 2等分，分点xi=a+（b-a）/ 2 ·i （i =0，1，2 · · · 2k）每个小区间长度为（b-a）/ 2 .2、参考Romberg数值积分，实现积分的数值求解，完成下列题目：三、算法设计及实现1、算法设计(a)function y=fun1(x)y=sqrt(4-(sin(x)).^2)；(b)function y=fun2 (x)y=log(1+2*x)/3;(c)function y=fun3(x)y=exp(-x).*sin(x.^2)四、实验步骤1、○1打开matlab软件，新建Romberg.m文件,在窗口中编辑Romber数值积分函数程序代码,并保存在指定的文件夹下,在Current Directory窗口右边点击《Browse For Folder》按钮指向Romberg.m文件；○2在Command Window中编辑相应要计算的题目的数值函数及相应的题目的表达式。
2、输出结果和初步分析说明（见附件一）。
五、使用说明实验结果分析1、在Command Window窗口中编辑要调用的函数名与指定的函数名字不同导致出现错误，通过改正与函数名相同即可。
2、Romberg方法也称为逐次分半加速法。
它是在梯形公式、辛卜生公式和柯特斯公式之间的关系的基础上，构造出一种加速计算积分的方法。
作为一种外推算法，它在不增加计算量的前提下提高了误差的精度在等距基点的情况下，用计算机计算积分值通常都采用把区间逐次分半的方法进行,这样，前一次分割得到的函数值在分半以后仍可被利用，且易于编程。
六、算法的改进和实验总结1、算法进一步的发展：在实际计算中，为防止假收敛，可设置最小二分次数k1，当k<k1时，跳过精度判别而继续运算。
2、用数值计算的方法来实现数值积分的求解，被调用的函数需要清楚的描述所要执行的问题的求解过程，在matlab函数调用时，执行函数一定要与函数同名。
七、源程序（见附件二）附件一>>[val,M]=Romberg('fun1',[0 pi/2],1e-10)val =2.9349M =2.9311 0 0 0 0 0 02.9349 2.9362 0 0 0 0 02.9349 2.9349 2.9348 0 0 0 02.9349 2.9349 2.9349 2.9349 0 0 02.9349 2.9349 2.9349 2.9349 2.9349 0 02.9349 2.9349 2.9349 2.9349 2.9349 2.9349 02.9349 2.9349 2.9349 2.9349 2.9349 2.9349 2.9349>>[val,M]=Romberg('fun2',[1 3],1e-10)y =0.3662 0.6486y =0.5365y =0.4621 0.5973y =0.4176 0.5014 0.5682 0.6239y =0.3929 0.4406 0.4823 0.5194 0.5527 0.5831 0.6109 0.6365y =0.3798 0.4055 0.4293 0.4515 0.4724 0.4920 0.5105 0.5280 0.5447 0.5606 0.5757 0.5902 0.6041 0.6175 0.6303 0.6426y =Columns 1 through 160.3731 0.3864 0.3992 0.4116 0.4235 0.4350 0.4461 0.4568 0.4673 0.4774 0.4872 0.4967 0.5060 0.5150 0.5237 0.5323Columns 17 through 320.5406 0.5488 0.5567 0.5644 0.5720 0.5794 0.5867 0.5938 0.6007 0.6075 0.6142 0.6207 0.6271 0.6334 0.6396 0.6456val =1.0543M =1.0148 0 0 0 0 0 01.0439 1.0536 0 0 0 0 01.0516 1.0542 1.0542 0 0 0 01.0536 1.0543 1.0543 1.0543 0 0 01.0541 1.0543 1.0543 1.0543 1.0543 0 01.0542 1.0543 1.0543 1.0543 1.0543 1.0543 01.0542 1.0543 1.0543 1.0543 1.0543 1.0543 1.0543>> [val,M]=Romberg('fun3',[0 2],1e-10)y = 0 -0.1024y = 0.3096y =0.1501 0.1736y =0.0486 0.2519 0.2865 0.0137y =0.0138 0.0963 0.2038 0.2889 0.3097 0.2400 0.0946 -0.0560y =0.0037 0.0291 0.0713 0.1228 0.1773 0.2289 0.2721 0.3016 0.3124 0.3011 0.2660 0.2089 0.1351 0.0536 -0.0234 -0.0828y =Columns 1 through 160.0009 0.0080 0.0209 0.0384 0.0596 0.0836 0.1094 0.1364 0.1637 0.1907 0.2166 0.2407 0.2624 0.2810 0.2958 0.3062 Columns 17 through 320.3117 0.3117 0.3061 0.2945 0.2770 0.2537 0.2251 0.1917 0.1547 0.1150 0.0740 0.0334 -0.0053 -0.0403 -0.0702 -0.0936y =Columns 1 through 160.0002 0.0021 0.0056 0.0107 0.0172 0.0249 0.0337 0.0434 0.0541 0.0654 0.0774 0.0899 0.1028 0.1161 0.1296 0.1432 Columns 17 through 320.1569 0.1705 0.1840 0.1973 0.2102 0.2228 0.2349 0.2464 0.2573 0.2674 0.2767 0.2851 0.2925 0.2988 0.3040 0.3080 Columns 33 through 480.3108 0.3122 0.3122 0.3109 0.3081 0.3038 0.2980 0.2907 0.2819 0.2717 0.2600 0.2470 0.2327 0.2171 0.2005 0.1828 Columns 49 through 640.1642 0.1449 0.1251 0.1048 0.0843 0.0638 0.0435 0.0235 0.0041 -0.0144 -0.0320 -0.0484 -0.0633 -0.0767 -0.0884 -0.0982val =0.2971M =-0.1024 0 0 0 0 0 0 00.2583 0.3786 0 0 0 0 0 00.2910 0.3019 0.2968 0 0 0 0 00.2957 0.2973 0.2970 0.2970 0 0 0 00.2967 0.2971 0.2971 0.2971 0.2971 0 0 00.2970 0.2971 0.2971 0.2971 0.2971 0.2971 0 00.2970 0.2971 0.2971 0.2971 0.2971 0.2971 0.2971 00.2971 0.2971 0.2971 0.2971 0.2971 0.2971 0.2971 0.2971 >>附件二：function [val,M]=Romberg(f,ab,dalt)% Romberg算法计算积分% [val,M]=Romberg(f,ab,dalt)% val 返回积分值% M T数表% f 被积函数,函数文件名% ab 积分区间% dalt 精度if nargin<3dalt=1e-7; %设置默认精度10-7endi=1; h=ab(2)-ab(1); t=0; j=0;T(1,1)=(h/2)*(feval(f,ab)*ones(2,1));while i<50 %最大跌代次数为50次i=i+1; h=h/2;T=[T zeros(i-1,1);zeros(1,i)];x=ab(1)+h:2*h:ab(2)-h;y=feval(f,x)*ones(size(x')); % feval(f,x)求得f在x点值T(i,1)=T(i-1,1)/2+h*y; % 细化区间,求得梯形值for t=2:ij=t-1;T(i,t)=(4^j*T(i,j)-T(i-1,j))/(4^j-1); %外推加速endif abs(T(i,i)-T(i-1,i-1))<=dalt %控制精度breakendendif nargout==2M=T; %按需求返回T数表endval=T(i,i); %积分值附件二（迭代法）：function val=DieDai(f,x0,delta)%用if nargin<3delta=1e-7;endn=1;x1=feval(f,x0);while abs(x1-x0)>deltax0=x1;x1=feval(f,x0);n=n+1;if n>100exit('跌代出问题');endendval=x1;disp(['跌代次数']);%显示跌代次数disp(n);附件二（牛顿法）：function val=NewtonDieDai(f,df,x0,delta)if nargin<4delta=1e-7;endn=1;x1=x0-feval(f,x0)/feval(df,x0);while abs(x1-x0)>deltax0=x1;x1=x0-feval(f,x0)/feval(df,x0)n=n+1if n>20exit('跌代可能不收敛！！！'); endendval=x1;disp(['跌代次数']);disp(n);。