cdq872 发表于 2007-7-14 10:37

灰色系统GM(1,1)-----灾变模型

CC ************************************
!灰色系统GM(1,1)-----灾变模型
cThe data is the rainfall of JINHUA from 1953-1986,and follow is the main program.
       program main
      use msimsl
      parameter(N=34,valve=700.0,M=7,Long=8)
      REAL x0(N),p(M),pp(M),p1(M),forcastp(Long+1),y(6,1),AINV(2,2),   
   & medium(2,M-1),b(M-1,2),a(2,M-1),u(2,1)
      data (x0(i),i=1,34)/336,1171,792,444,311,663,553,501,554,710,
   &                  509,424,536,437,720,645,723,607,573,384,
   &      917,533,535,627,773,333,360,402,344,439,
   &                  631,660,411,392/
c   Take the catastrophic sequence   
         k=1
       do i=1,N
      if(x0(i)>=valve)then
         p(k)=i   !时间序列
         k=k+1
      end if
       end do
WRITE(*,*) 'P=',P
c   Take the first accumulation and make the matrix B and Y.
          p1(1)=p(1)
       do i=2,M
          p1(i)=p1(i-1)+p(i) !把>=700数序号累加生成时间新序列
       end do
       do i=2,M
          y(i-1,1)=p(i)
       end do
       do i=1,M-1
      b(i,1)=-(p1(i)+p1(i+1))/2!B矩阵,紧邻均值生成B(M-1,2)
      b(i,2)=1
       end do
C   --------------------   
          A=TRANSPOSE(b)    !B矩阵转置:A(2,M-1)
AINV=MATMUL(A,B)   !矩阵乘法:AINV=BT*B: AINV(2,2)
WRITE(*,*) 'AINV=',AINV
    CALL INVA(AINV,2)!矩阵求逆   :AINV-1(2,2)
WRITE(*,*) 'AINV=',AINV      
!medium=(BT*B)-1*BT   :medium(2,M-1)
      medium=MATMUL(AINV,TRANSPOSE(b)) !TRANSPOSE矩阵转置MATMUL矩阵乘法
!u=(BT*B)-1*BT*y      :u(2,1)
       u=MATMUL(medium,y)!最小二乘法求灰参数
C   --------------------   
       forcastp(1)=p1(1)
         do i=1,Long   !计算微分方程的解,i+1为预测值
       forcastp(i+1)=(p1(1)-u(2,1)/u(1,1))*exp(-u(1,1)*i)+u(2,1)/u(1,1)
         end do
         pp(1)=P1(1)
       do i=2,M
      pp(i)=forcastp(i)-forcastp(i-1)   !原始序列的预测值
       end do
       write(*,*)'forcastp=',forcastp
       write(*,*)'原序列预测值pp＝',pp
       write(*,*)'The next catastrophic year is:'
       write(*,*) 1952+(forcastp(Long)-forcastp(Long-1))
       pause
      end
C___________INVA()：求逆 ______A-1(2,2)_______
      SUBROUTINE INVA(A,M)
      DIMENSION A(M,M)
C       ---------------
      DO K=1,M
          DKK=A(K,K)
       A(K,K)=1.0
         DO J=1,M
       A(K,J)=A(K,J)/DKK
         ENDDO
       DO 13190 I=1,M
      IF (I.EQ.K) GO TO 13190
       DIK=A(I,K)
       A(I,K)=0
          DO J=1,M
       A(I,J)=A(I,J)-DIK*A(K,J)
          ENDDO
13190   CONTINUE
      ENDDO
      RETURN
      END

cdq872 发表于 2007-7-14 16:56

以上为灾变模型，程序在cvf6.6通过。
尚有灰色系统预测模型，有机会上传。


求rbf的fortran程序，请"风花雪月"版主指点。

也请路过的高人指点，
邮箱cdq872@163.com，或加我QQ：136715186

风花雪月 发表于 2007-7-23 15:50

原帖由 cdq872 于 2007-7-14 16:56 发表 http://www.chinavib.com/forum/images/common/back.gif
以上为灾变模型，程序在cvf6.6通过。
尚有灰色系统预测模型，有机会上传。

希望能够尽快看到

风花雪月 发表于 2007-7-23 15:52

原帖由 cdq872 于 2007-7-14 16:56 发表 http://www.chinavib.com/forum/images/common/back.gif
求rbf的fortran程序，请"风花雪月"版主指点。

也请路过的高人指点，
邮箱cdq872@163.com，或加我QQ：136715186

http://orion.math.iastate.edu/burkardt/f_src/rbf/rbf.html

cdq872 发表于 2007-7-25 15:08

谢谢"风花雪月"
上面地址中的程序缺关键子程序,所以看不懂.

风花雪月 发表于 2007-7-29 15:54

原帖由 cdq872 于 2007-7-25 15:08 发表 http://www.chinavib.com/forum/images/common/back.gif
谢谢"风花雪月"
上面地址中的程序缺关键子程序,所以看不懂.

注意看页面说明

RBF references the MINPACK library.
页: [1]
查看完整版本: 灰色系统GM(1,1)-----灾变模型