马上注册，结交更多好友，享用更多功能，让你轻松玩转社区。
您需要 登录 才可以下载或查看，没有账号？我要加入
x
CC ************************************
! 灰色系统GM(1,1)-----灾变模型
c The 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 |
发表于 2007-7-14 10:37
提升卡
变色卡
显身卡
楼主
