Romberg数值积分.docx

Romberg数值积分实验报告一、实验目的1、加深外推法的原理理解, 掌握Romberg外推法的计算方法。2、用matlab软件实现Romberg数值积分来计算题目的运算。二、基本理论及背景1、理论推导：对区间a， b，令h=b-a构造梯形值序列T2K。 T1=hf(a)+f(b)/2 把区间二等分，每个小区间长度为 h/2=(b-a)/2，于是 T2 =T1/2+h/2f(a+h/2） 把区间四(2)等分，每个小区间长度为h/2 =（b-a）/4，于是 T4 =T2/2+h/2f（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、打开matlab软件，新建Romberg.m文件,在窗口中编辑Romber数值积分函数程序代码,并保存在指定的文件夹下,在Current Directory窗口右边点击Browse For Folder按钮指向Romberg.m文件；在Command Window中编辑相应要计算的题目的数值函数及相应的题目的表达式。2、输出结果和初步分析说明（见附件一）。五、使用说明实验结果分析1、在Command Window窗口中编辑要调用的函数名与指定的函数名字不同导致出现错误，通过改正与函数名相同即可。2、Romberg方法也称为逐次分半加速法。它是在梯形公式、辛卜生公式和柯特斯公式之间的关系的基础上，构造出一种加速计算积分的方法。作为一种外推算法，它在不增加计算量的前提下提高了误差的精度在等距基点的情况下，用计算机计算积分值通常都采用把区间逐次分半的方法进行,这样，前一次分割得到的函数值在分半以后仍可被利用，且易于编程。六、算法的改进和实验总结 1、算法进一步的发展：在实际计算中，为防止假收敛，可设置最小二分次数k1，当k val,M=Romberg(fun1,0 pi/2,1e-10)val = 2.9349M = 2.9311 0 0 0 0 0 0 2.9349 2.9362 0 0 0 0 0 2.9349 2.9349 2.9348 0 0 0 0 2.9349 2.9349 2.9349 2.9349 0 0 0 2.9349 2.9349 2.9349 2.9349 2.9349 0 0 2.9349 2.9349 2.9349 2.9349 2.9349 2.9349 0 2.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 16 0.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.5323 Columns 17 through 32 0.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 0 1.0439 1.0536 0 0 0 0 0 1.0516 1.0542 1.0542 0 0 0 0 1.0536 1.0543 1.0543 1.0543 0 0 0 1.0541 1.0543 1.0543 1.0543 1.0543 0 0 1.0542 1.0543 1.0543 1.0543 1.0543 1.0543 0 1.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 16 0.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 32 0.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 16 0.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 32 0.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 48 0.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 64 0.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 0 0.2583 0.3786 0 0 0 0 0 0 0.2910 0.3019 0.2968 0 0 0 0 0 0.2957 0.2973 0.2970 0.2970 0 0 0 0 0.2967 0.2971 0.2971 0.2971 0.2971 0 0 0 0.2970 0.2971 0.2971 0.2971 0.2971 0.2971 0 0 0.2970 0.2971 0.2971 0.2971 0.2971 0.2971 0.2971 0 0.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 nargin3 dalt=1e-7; %设置默认精度10-7 endi=1; h=ab(2)-ab(1); t=0; j=0;T(1,1)=(h/2)*(feval(f,ab)*ones(2,1);while i50 %最大跌代次数为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:i j=t-1; T(i,t)=(4j*T(i,j)-T(i-1,j)/(4j-1); %外推加速 end if abs(T(i,i)-T(i-1,i-1)=dalt %控制精度 break endendif nargout=2 M=T; %按需求返回T数表endval=T(i,i); %积分值附件二（迭代法）：function val=DieDai(f,x0,delta)%用if nargindelta x0=x1; x1=feval(f,x0); n=n+1; if n100