romberg1.m

来自「高频地波雷达完整仿真matlab程序 网上找到的」· M 代码 · 共 27 行

function [I,R]=Romberg1(a,b,eps2)
%f表示被积函数；
%a,b表示被积区间；
%eps2表示控制精度；
%I表示用Romberg加速算法所求得的积分值；
%R为Romberg表。
 h=b-a;
 R(1,1)=h*(f(a)+f(b))/2;
 M=1;
 J=0;
 err=1;
 while err>eps2
     J=J+1;
     h=h/2;
     S=0;
     for p=1:M
         x=a+h*(2*p-1);
         S=S+f(x);
     end
     R(J+1,1)=R(J,1)/2+h*S;
     M=M*2;
     for K=1:J
         R(J+1,K+1)=R(J+1,K)+(R(J+1,K)-R(J,K))/(4^K-1);
     end
     err=abs(R(J+1,J)-R(J+1,J+1));
 end
 I=R(J+1,J+1);