Demo entry 6645807

py

   

Submitted by 张樱潇 on Oct 11, 2017 at 17:50
Language: Fortran. Code size: 4.1 kB.

program zyx04
implicit none
    real,dimension(10)::x1,x2,x3,x4,y,x0,y0    !x0,y0是原始自变量,x1,x2,x3,x4,y为变换后的自变量,转换为多元线性拟合
    integer,parameter::n=4                     !n为变量的个数
    real t,t1,t2,t3,t4                         !分别为各项代换后变量的平均值
    real b0                                    !常系数
    real::s(n,n)=0.0                           !xi,xj之间的协方差矩阵
    real,dimension(n)::ss=0.0,b                !ss为xi,y之间的协方差矩阵,b为系数矩阵
    integer i,j
    real::w=0.0                                   !求出拟合函数与原始点之间的误差

    x0(1:10)=(/1.02,0.95,0.87,0.77,0.67,0.56,0.44,0.30,0.16,0.01/)
    y0(1:10)=(/0.39,0.32,0.27,0.22,0.18,0.15,0.13,0.12,0.13,0.15/)

    do i=1,10                                  !为每个新定义变量的数组赋值
        x1(i)=x0(i)
        x2(i)=y0(i)
        x3(i)=x0(i)*y0(i)
        x4(i)=y0(i)*y0(i)
        y(i)=x0(i)*x0(i)
    end do

    t=0
    t1=0
    t2=0
    t3=0
    t4=0

    do i=1,10                                   !求新定义变量的平均值
        t1=t1+x1(i)
        t2=t2+x2(i)
        t3=t3+x3(i)
        t4=t4+x4(i)
        t=t+y(i)
    end do

    t1=t1/10
    t2=t2/10
    t3=t3/10
    t4=t4/10
    t=t/10

    do j=1,10                                   !求xi,xj之间的协方差
        s(1,1)=s(1,1)+(x1(j)-t1)*(x1(j)-t1)
        s(1,2)=s(1,2)+(x1(j)-t1)*(x2(j)-t2)
        s(1,3)=s(1,3)+(x1(j)-t1)*(x3(j)-t3)
        s(1,4)=s(1,4)+(x1(j)-t1)*(x4(j)-t4)
        s(2,2)=s(2,2)+(x2(j)-t2)*(x2(j)-t2)
        s(2,3)=s(2,3)+(x2(j)-t2)*(x3(j)-t3)
        s(2,4)=s(2,4)+(x2(j)-t2)*(x4(j)-t4)
        s(3,3)=s(3,3)+(x3(j)-t3)*(x3(j)-t3)
        s(3,4)=s(3,4)+(x3(j)-t3)*(x4(j)-t4)
        s(4,4)=s(4,4)+(x4(j)-t4)*(x4(j)-t4)
    end do
    s(2,1)=s(1,2)
    s(3,1)=s(1,3)
    s(4,1)=s(1,4)
    s(3,2)=s(2,3)
    s(4,2)=s(2,4)
    s(4,3)=s(3,4)

    do i=1,4
        do j=1,4
            s(i,j)=s(i,j)/10
        end do
    end do

    do i=1,10                                  !求xi,y之间的协方差
        ss(1)=ss(1)+(x1(i)-t1)*(y(i)-t)
        ss(2)=ss(2)+(x2(i)-t2)*(y(i)-t)
        ss(3)=ss(3)+(x3(i)-t3)*(y(i)-t)
        ss(4)=ss(4)+(x4(i)-t4)*(y(i)-t)
    end do

    do i=1,4
        ss(i)=ss(i)/10
    end do

    write(*,*)"协方差矩阵"
    do i=1,4
        write(*,*)s(i,:)
    end do

    write(*,*)"协方差逆矩阵"
    call ni(s,n)                              !调用函数,求s的逆矩阵
    do i=1,4
        write(*,*)s(i,:)
    end do

    b=matmul(s,ss)
    b0=t-t1*b(1)-t2*b(2)-t3*b(3)-t4*b(4)

    write(*,*)"系数"
    write(*,*)b0,b

    do i=1,10
        w=w+(b0+b(1)*x1(i)+b(2)*x2(i)+b(3)*x3(i)+b(4)*x4(i)-y(i))**2
    end do

    write(*,*)"误差"
    write(*,*)w

    open(1,file='4.txt',form='formatted',status='replace')
    write(1,*)"系数"
    write(1,*)b0,b
    write(1,*)"误差"
    write(1,*)w
    close(1)

stop
end program

subroutine ni(A,N)    !求解逆矩阵的子程序,子程序参考徐良士算法集brinv逆矩阵求解算法
implicit none
integer::N
integer::I,J,K
real::A(N,N),IS(N),JS(N)
real::T,D,L
	L=1
	DO 100 K=1,N
	  D=0.0
	  DO 10 I=K,N
	  DO 10 J=K,N
	    IF (ABS(A(I,J)).GT.D) THEN
	      D=ABS(A(I,J))
	      IS(K)=I
	      JS(K)=J
	    END IF
10	  CONTINUE
	  IF (D+1.0.EQ.1.0) THEN
	    L=0
	    WRITE(*,20)
	    RETURN
	  END IF
20	  FORMAT(1X,'ERR**NOT INV')
	  DO 30 J=1,N
	    T=A(K,J)
	    A(K,J)=A(IS(K),J)
	    A(IS(K),J)=T
30	  CONTINUE
	  DO 40 I=1,N
	    T=A(I,K)
	    A(I,K)=A(I,JS(K))
	    A(I,JS(K))=T
40	  CONTINUE
	  A(K,K)=1/A(K,K)
	  DO 50 J=1,N
	    IF (J.NE.K) THEN
	      A(K,J)=A(K,J)*A(K,K)
	    END IF
50	  CONTINUE
	  DO 70 I=1,N
	    IF (I.NE.K) THEN
	      DO 60 J=1,N
	        IF (J.NE.K) THEN
	          A(I,J)=A(I,J)-A(I,K)*A(K,J)
	        END IF
60	      CONTINUE
	    END IF
70	  CONTINUE
	  DO 80 I=1,N
	    IF (I.NE.K) THEN
	      A(I,K)=-A(I,K)*A(K,K)
	    END IF
80	  CONTINUE
100	CONTINUE
	DO 130 K=N,1,-1
	  DO 110 J=1,N
	    T=A(K,J)
	    A(K,J)=A(JS(K),J)
	    A(JS(K),J)=T
110	  CONTINUE
	  DO 120 I=1,N
	    T=A(I,K)
	    A(I,K)=A(I,IS(K))
	    A(I,IS(K))=T
120	  CONTINUE
130	CONTINUE
	RETURN
	END

This snippet took 0.01 seconds to highlight.

Back to the Entry List or Home.

Delete this entry (admin only).