rombergfunc.m

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

function [A]=Romberg(sitaL,sitaR,eta,u,q,f0)
%Romberg积分法
%sitaL,sitaR表示被积区间；
%A表示用Romberg加速算法所求得的积分值；
%eta,u,q,f0风别为归一化多普勒频率，风速，风向及雷达工作频率。
 eps2=10^-7;
 h=sitaR-sitaL;
 R(1,1)=h*(integrandfunc(eta,sitaL,u,f0,q)+integrandfunc(eta,sitaR,u,f0,q))/2;
 M=1;
 J=0;
 err=1;
 while err>eps2
     J=J+1;
     h=h/2;
     S=0;
     for p=1:M
         sita=sitaL+h*(2*p-1);
         S=S+integrandfunc(eta,sita,u,f0,q);
     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
 A=R(J+1,J+1);