pts.m

自适应梯形求积法 特别适合被积函数为强峰的定积分 此程序修改余地很大

% 自适应梯形求积法  特别适合被积函数为强峰的定积分  此程序修改余地很大
% 1.积分区间 [a,b] 分成两个相等的区间，然后，每个子区间在分
% 2.判断  | t01-（t02+t12）| < eps/1.4  | t11-（t22+t32）| < eps/1.4  近似积分值为 t ＝ t02+t12+t22+t32
% 3.所有子区间都要满足精度要求 



% SUBROUTINE PTS(A,B,F,EPS,T)
% A   实型变量，输入参数。积分下限
% B   实型变量，输入参数。积分上限
% F   实型函数段名，输入参数。 用于计算被积函数的函数值。与之对应的实函数段名在主程序中用外部语句说明。此函数子程序由用户自编。其语句形式为 FUNCTION F（X）
% T   实型变量，输出参数。存放积分值。
% EPS 实型变量，输入参数。精确度。

    SUBROUTINE SIMP(A,B,F,EPS,T)
    DIMENSION S(30,7)
    DOUBLE PRECISION A,B,F,T,S,F0,F1,P,H,X,F3,T1,T2
 
    T=0.0
    F0=F(A)
    F1=F(B)
    K=1
    S(K,1)=A
    S(K,2)=B
    S(K,3)=B-A
    S(K,4)=F0
    S(K,5)=F1
    S(K,6)=P
    S(K,7)=EPS



10   IF (K.NE.0) THEN
         H=S(K,3)
         X=S(K,1)+H/2.0
         F3=F(X)
         T1=H*(S(K,4)+F3)/4.0
         T2=H*(S(K,5)+F3)/4.0
      
          IF ((ABS(S(K,6)-T1-T2).LT.S(K,7).OR.(K.GE.29))THEN
               T=T+T1+T2
               K=K-1
          ELSE

               S(K+1,1)=X
               S(K+1,2)=S(K,2)
               S(K,2)=X
               S(K,3)=H/2.0
               S(K+1,3)=S(K,3)
               S(K+1,4)=F3
               S(K+1,5)=S(K,5)
               S(K,5)=F3
               S(K,7)=S(K,7)/1.4
               S(K+1,7)=S(K,7)
               S(K,6)=T1
               S(K+1,6)=T2
               K=K+1
           END IF
           GOTO 10
      END IF
      
      RETURN

      END