fortran 程序(应该是很简单的小程序)

fortran 程序(应该是很简单的小程序),第1张

Fortran源自于“公式翻译”(英语:FormulaTranslation)的缩写,是一种编程语言

它是世界上芦肢世最早出现的计算机高级程序设计语言,广泛应用于科学和工程计算领域。FORTRAN语言以其特有的功能在数值、科学和工程计算领域发挥着重要作用。

随着FORTRAN语言版本的不断更新和变化,饥旅语言不兼容性问题日益突出,语言标准化工作被提上了日程。

1962年5月:美国标准化协会(简称ANSI)着手进行FORTRAN语言标准化的研究工作。

1966年:ANSI正式公布了两个标准文本:美国国家标准FORTRAN(ANSI X3.9-1966)和美国国家标准基本FORTRAN(ANSI X3.10-1966),前者相当于FORTRAN Ⅳ,后者相当于FORTRANⅡ。基本FORTRAN是美国国家标准FORTRAN的一个子集,从而实现了语言的向下兼容,初步解决了语言的兼陪肢容性问题。

!高纤册斯消去法

subroutine agaus(a,b,n,x,l,js)

dimension a(n,n),x(n),b(n),js(n)

double precision a,b,x,t

l=1 !逻辑变量

do k=1,n-1

d=0.0

do i=k,n

do j=k,n

if (abs(a(i,j))>d) then

d=abs(a(i,j))

js(k)=j

is=i

end if

end do

end do !把行绝对值最大的元素换到粗蠢主元位毁凳宏置

if (d+1.0==1.0) then

l=0

else !最大元素为0无解

if(js(k)/=k) then

do i=1,n

t=a(i,k)

a(i,k)=a(i,js(k))

a(i,js(k))=t

end do !最大元素不在K行,K行

end if

if(is/=k) then

do j=k,n

t=a(k,j)

a(k,j)=a(is,j)

a(is,j)=t !交换到K列

end do

t=b(k)

b(k)=b(is)

b(is)=t

end if !最大元素在主对角线上

end if !消去

if (l==0) then

write(*,100)

return

end if

do j=k+1,n

a(k,j)=a(k,j)/a(k,k)

end do

b(k)=b(k)/a(k,k) !求三角矩阵

do i=k+1,n

do j=k+1,n

a(i,j)=a(i,j)-a(i,k)*a(k,j)

end do

b(i)=b(i)-a(i,k)*b(k)

end do

end do

if (abs(a(n,n))+1.0==1.0) then

l=0

write(*,100)

return

end if

x(n)=b(n)/a(n,n)

do i=n-1,1,-1

t=0.0

do j=i+1,n

t=t+a(i,j)*x(j)

end do

x(i)=b(i)-t

end do

100 format(1x,'fail')

js(n)=n

do k=n,1,-1

if (js(k)/=k) then

t=x(k)

x(k)=x(js(k))

x(js(k))=t

end if

end do

return

end

program main

dimension a(4,4),b(4),x(4),js(4)

double precision a,b,x

real m1,m2,j

open(1,file="laiyi.txt")

read(1,*)m1,m2,j

close(1)

n=4

print*,m1,m2,j

a(1,1)=m1*cos(3.14159*j/180)

a(1,2)=-m1

a(1,3)=-sin(3.14159*j/180)

a(1,4)=0

a(2,1)=m1*sin(3.14159*j/180)

a(2,2)=0

a(2,3)=cos(3.14159*j/180)

a(2,4)=0

a(3,1)=0

a(3,2)=m2

a(3,3)=-sin(3.14159*j/180)

a(3,4)=0

a(4,1)=0

a(4,2)=0

a(4,3)=-cos(3.14159*j/180)

a(4,4)=1

b(1)=0

b(2)=m1*9.8

b(3)=0

b(4)=m2*9.8

call agaus(a,b,n,x,l,js)

if (l/=0) then

write(*,*)"a1=",x(1),"a2=",x(2) ,"n1=",x(3),"n2=",x(4)

end if

end

!逆矩阵求解

SUBROUTINE qiuni(A,N,L,IS,JS)

DIMENSION A(N,N),IS(N),JS(N)

DOUBLE PRECISION A,T,D

L=1

DO K=1,N

D=0.0

DO I=K,N

DO J=K,N

IF(ABS(A(I,J)).GT.D) THEN !把最大的元素给D

D=ABS(A(I,J))

IS(K)=I

JS(K)=J

END IF

END DO

END DO

IF (D+1.0.EQ.1.0)THEN

L=0

WRITE(*,200)

RETURN

END IF

200 FORMAT(1X,'ERR**NOT INV')

DO J=1,N

T=A(K,J)

A(K,J)=A(IS(K),J)

A(IS(K),J)=T

END DO

DO I=1,N

T=A(I,K)

A(I,K)=A(I,JS(K))

A(I,JS(K))=T

END DO

A(K,K)=1/A(K,K)

DO J=1,N

IF(J.NE.K)THEN

A(K,J)=A(K,J)*A(K,K)

END IF

END DO

DO I=1,N

IF(I.NE.K)THEN

DO J=1,N

IF(J.NE.K)THEN

A(I,J)=A(I,J)-A(I,K)*A(K,J)

END IF

END DO

END IF

END DO

DO I=1,N

IF(I.NE.K)THEN

A(I,K)=-A(I,K)*A(K,K)

END IF

END DO

END DO

DO K=N,1,-1

DO J=1,N

T=A(K,J)

A(K,J)=A(JS(K),J)

A(JS(K),J)=T

END DO

DO I=1,N

T=A(I,K)

A(I,K)=A(I,IS(K))

A(I,IS(K))=T

END DO

END DO

RETURN

END

SUBROUTINE BRMUL(A,B,N,C)

DIMENSION A(N,N),B(N),C(N)

DOUBLE PRECISION A,B,C

DO I=1,N

DO J=1,N

C(I)=0.0

DO L=1,N

C(I)=C(I)+A(I,L)*B(L)

END DO

END DO

END DO

RETURN

END

program main

DIMENSION A(4,4),B(4,1),C(4,1),IS(4),JS(4)

DOUBLE PRECISION A,B,C

REAL M1,M2,JD

OPEN(1,FILE='LAIYI.TXT')

READ(1,*) M1,M2,JD

PRINT*,M1,M2,JD

CLOSE(1)

A(1,1)=M1*COS(3.14*JD/180)

A(1,2)=-M1

A(1,3)=-SIN(3.14*JD/180)

A(1,4)=0

A(2,1)=M1*SIN(3.14*JD/180)

A(2,2)=0

A(2,3)=COS(3.14*JD/180)

A(2,4)=0

A(3,1)=0

A(3,2)=M2

A(3,3)=-SIN(3.14*JD/180)

A(3,4)=0

A(4,1)=0

A(4,2)=0

A(4,3)=-COS(3.14*JD/180)

A(4,4)=1

B(1,1)=0

B(2,1)=M1*9.8

B(3,1)=0

B(4,1)=M2*9.8

CALL QIUNI(A,4,L,IS,JS)

CALL BRMUL(A,B,4,C)

WRITE(*,*) (C(I,1),I=1,4)

END

画图

USE MSFLIB

INTEGER status

TYPE(xycoord) xy

status=SETCOLORRGB(#FFFFFF)

status1=SETCOLORRGB(#0000FF)

OPEN(1,FILE="G.TXT")

READ(1,*) G1,G2,G3,G4

OPEN(2,FILE="N.TXT")

READ(2,*) N1,N2,N3,N4

CALL MOVETO(INT(20),INT(20),xy)

status=LINETO(INT(40),INT(G1))

status=LINETO(INT(80),INT(G2))

status=LINETO(INT(120),INT(G3))

status=LINETO(INT(160),INT(G4))

CALL SETLINESTYLE(#FF00)

CALL MOVETO(INT(20),INT(20),xy)

status1=LINETO(INT(40),INT(N1))

status1=LINETO(INT(80),INT(N2))

status1=LINETO(INT(120),INT(N3))

status1=LINETO(INT(160),INT(N4))

READ(*,*)

END

如果对您有帮助,请记得采纳为满意答案,谢谢!祝您生活愉快!

module data_type

implicit none

integer(kind=4), parameter :: IB=4, RP=8

end module data_type

module data_Rosen

use data_type

implicit none

integer(kind=IB), parameter :: Dim_XC=10

end module data_Rosen

module data_HDE

use data_type

use data_Rosen

implicit none

integer(kind=IB), parameter :: NP=20, itermax=20000, strategy=6, &

refresh=500, iwrite=7

integer(kind=IB), dimension(3), parameter :: method=(/0, 1, 0/)

real(kind=RP), parameter :: VTR=-1.0e-4_RP, CR_XC=0.5_RP

real(kind=RP) :: F_XC=0.8_RP, F_CR=0.8_RP

real(kind=RP), dimension(Dim_XC), parameter :: XCmin=-10.0_RP, XCmax=10.0_RP

real(kind=RP), dimension(Dim_XC) :: bestmem_XC

integer(kind=IB) :: nfeval

real(kind=RP) :: bestval

end module data_HDE

program Rosen

use data_type

use data_Rosen

use data_HDE

implicit none

integer(kind=IB) :: i

integer (kind=IB), dimension(8) :: time

intrinsic date_and_time

external FTN

open(iwrite,file='Rosen.txt')

call date_and_time(values=time)

write(unit=iwrite, FMT=11) time(1:3), time(5:7)

call DE_Fortran90(FTN, Dim_XC, XCmin, XCmax, VTR, NP, itermax, F_XC,&握毁

CR_XC, strategy, refresh, iwrite, bestmem_XC, &

bestval, nfeval, F_CR, method)

write(iwrite,205) NP, nfeval, method(1:3)

write(iwrite,FMT=201) F_XC, CR_XC, F_CR

write(iwrite,FMT=200) bestval

do i=1,Dim_XC

write(iwrite,FMT=202) i,bestmem_XC(i)

end do

200 format(/2x, 'Bestval=', ES14.7)

201 format(2x, 'F_XC =',F6.3, 2x, 'CR_XC =', F6.3, 2x, '尘桥F_CR =', F6.3)

202 format(2x, 'best_XC(',I3,') =',ES14.7)

205 format(2x, 'NP=', I4, 4x, 'No. function call =', I9, &

/2x, 'mehtod(1:3) =',3I3)

call date_and_time(values=time)

write(unit=iwrite, FMT=10)time(1:3), time(5:7)

10 format(/1x, 'End of Program. Date:', I4, '/', I2,'/', I2, '段兄备, Time: ', I2,':',I2,':',I2)

11 format(1x, 'Beginning of Program. Date:', I4, '/', I2,'/', I2, ', Time: ', I2,':',I2,':',I2)

end program Rosen

subroutine FTN(X, objval)

use data_type

use data_Rosen

implicit none

real(kind=RP), dimension(Dim_XC), intent(in) :: X

real(kind=RP), intent(out) :: objval

integer(kind=IB) :: i

i=Dim_XC

objval=sum(100.0*(x(1:i-1)**2-x(2:i))**2+(1.0-x(1:i-1))**2)

end subroutine FTN

subroutine DE_Fortran90(obj, Dim_XC, XCmin, XCmax, VTR, NP, itermax, F_XC, &

CR_XC, strategy, refresh, iwrite, bestmem_XC, bestval, nfeval, &

F_CR, method)

!.......................................................................

!

! Differential Evolution for Optimal Control Problems

!

!.......................................................................

! This Fortran 90 program translates from the original MATLAB

! version of differential evolution (DE). This FORTRAN 90 code

! has been tested on Compaq Visual Fortran v6.1.

! Any users new to the DE are encouraged to read the article of Storn and Price.

!

! Refences:

! Storn, R., and Price, K.V., (1996). Minimizing the real function of the

! ICEC'96 contest by differential evolution. IEEE conf. on Evolutionary

! Comutation, 842-844.

!

! This Fortran 90 program written by Dr. Feng-Sheng Wang

! Department of Chemical Engineering, National Chung Cheng University,

! Chia-Yi 621, Taiwan, e-mail: chmfsw@ccunix.ccu.edu.tw

!.........................................................................

! obj : The user provided file for evlauting the objective function.

! subroutine obj(xc,fitness)

! where "xc" is the real decision parameter vector.(input)

! "fitness" is the fitness value.(output)

! Dim_XC : Dimension of the real decision parameters.

! XCmin(Dim_XC) : The lower bound of the real decision parameters.

! XCmax(Dim_XC) : The upper bound of the real decision parameters.

! VTR : The expected fitness value to reach.

! NP : Population size.

! itermax : The maximum number of iteration.

! F_XC : Mutation scaling factor for real decision parameters.

! CR_XC : Crossover factor for real decision parameters.

! strategy : The strategy of the mutation operations is used in HDE.

! refresh : The intermediate output will be produced after "refresh"

! iterations. No intermediate output will be produced if

! "refresh <1".

! iwrite : The unit specfier for writing to an external data file.

! bestmen_XC(Dim_XC) : The best real decision parameters.

! bestval : The best objective function.

! nfeval : The number of function call.

! method(1) = 0, Fixed mutation scaling factors (F_XC)

! = 1, Random mutation scaling factors F_XC=[0, 1]

! = 2, Random mutation scaling factors F_XC=[-1, 1]

! method(2) = 1, Random combined factor (F_CR) used for strategy = 6

! in the mutation operation

! = other, fixed combined factor provided by the user

! method(3) = 1, Saving results in a data file.

! = other, displaying results only.

use data_type, only : IB, RP

implicit none

integer(kind=IB), intent(in) :: NP, Dim_XC, itermax, strategy, &

iwrite, refresh

real(kind=RP), intent(in) :: VTR, CR_XC

real(kind=RP) :: F_XC, F_CR

real(kind=RP), dimension(Dim_XC), intent(in) :: XCmin, XCmax

real(kind=RP), dimension(Dim_XC), intent(inout) :: bestmem_XC

real(kind=RP), intent(out) :: bestval

integer(kind=IB), intent(out) :: nfeval

real(kind=RP), dimension(NP,Dim_XC) :: pop_XC, bm_XC, mui_XC, mpo_XC, &

popold_XC, rand_XC, ui_XC

integer(kind=IB) :: i, ibest, iter

integer(kind=IB), dimension(NP) :: rot, a1, a2, a3, a4, a5, rt

integer(kind=IB), dimension(4) :: ind

real(kind=RP) :: tempval

real(kind=RP), dimension(NP) :: val

real(kind=RP), dimension(Dim_XC) :: bestmemit_XC

real(kind=RP), dimension(Dim_XC) :: rand_C1

integer(kind=IB), dimension(3), intent(in) :: method

external obj

intrinsic max, min, random_number, mod, abs, any, all, maxloc

interface

function randperm(num)

use data_type, only : IB

implicit none

integer(kind=IB), intent(in) :: num

integer(kind=IB), dimension(num) :: randperm

end function randperm

end interface

!!-----Initialize a population --------------------------------------------!!

pop_XC=0.0_RP

do i=1,NP

call random_number(rand_C1)

pop_XC(i,:)=XCmin+rand_C1*(XCmax-XCmin)

end do

!!--------------------------------------------------------------------------!!

!!------Evaluate fitness functions and find the best member-----------------!!

val=0.0_RP

nfeval=0

ibest=1

call obj(pop_XC(1,:), val(1))

bestval=val(1)

nfeval=nfeval+1

do i=2,NP

call obj(pop_XC(i,:), val(i))

nfeval=nfeval+1

if (val(i) <bestval) then

ibest=i

bestval=val(i)

end if

end do

bestmemit_XC=pop_XC(ibest,:)

bestmem_XC=bestmemit_XC

!!--------------------------------------------------------------------------!!

bm_XC=0.0_RP

rot=(/(i,i=0,NP-1)/)

iter=1

!!------Perform evolutionary computation------------------------------------!!

do while (iter <= itermax)

popold_XC=pop_XC

!!------Mutation operation--------------------------------------------------!!

ind=randperm(4)

a1=randperm(NP)

rt=mod(rot+ind(1),NP)

a2=a1(rt+1)

rt=mod(rot+ind(2),NP)

a3=a2(rt+1)

rt=mod(rot+ind(3),NP)

a4=a3(rt+1)

rt=mod(rot+ind(4),NP)

a5=a4(rt+1)

bm_XC=spread(bestmemit_XC, DIM=1, NCOPIES=NP)

!----- Generating a random sacling factor--------------------------------!

select case (method(1))

case (1)

call random_number(F_XC)

case(2)

call random_number(F_XC)

F_XC=2.0_RP*F_XC-1.0_RP

end select

!---- select a mutation strategy-----------------------------------------!

select case (strategy)

case (1)

ui_XC=bm_XC+F_XC*(popold_XC(a1,:)-popold_XC(a2,:))

case default

ui_XC=popold_XC(a3,:)+F_XC*(popold_XC(a1,:)-popold_XC(a2,:))

case (3)

ui_XC=popold_XC+F_XC*(bm_XC-popold_XC+popold_XC(a1,:)-popold_XC(a2,:))

case (4)

ui_XC=bm_XC+F_XC*(popold_XC(a1,:)-popold_XC(a2,:)+popold_XC(a3,:)-popold_XC(a4,:))

case (5)

ui_XC=popold_XC(a5,:)+F_XC*(popold_XC(a1,:)-popold_XC(a2,:)+popold_XC(a3,:) &

-popold_XC(a4,:))

case (6) ! A linear crossover combination of bm_XC and popold_XC

if (method(2) == 1) call random_number(F_CR)

ui_XC=popold_XC+F_CR*(bm_XC-popold_XC)+F_XC*(popold_XC(a1,:)-popold_XC(a2,:))

end select

!!--------------------------------------------------------------------------!!

!!------Crossover operation-------------------------------------------------!!

call random_number(rand_XC)

mui_XC=0.0_RP

mpo_XC=0.0_RP

where (rand_XC <CR_XC)

mui_XC=1.0_RP

! mpo_XC=0.0_RP

elsewhere

! mui_XC=0.0_RP

mpo_XC=1.0_RP

end where

ui_XC=popold_XC*mpo_XC+ui_XC*mui_XC

!!--------------------------------------------------------------------------!!

!!------Evaluate fitness functions and find the best member-----------------!!

do i=1,NP

!!------Confine each of feasible individuals in the lower-upper bound-------!!

ui_XC(i,:)=max(min(ui_XC(i,:),XCmax),XCmin)

call obj(ui_XC(i,:), tempval)

nfeval=nfeval+1

if (tempval <val(i)) then

pop_XC(i,:)=ui_XC(i,:)

val(i)=tempval

if (tempval <bestval) then

bestval=tempval

bestmem_XC=ui_XC(i,:)

end if

end if

end do

bestmemit_XC=bestmem_XC

if( (refresh >0) .and. (mod(iter,refresh)==0)) then

if (method(3)==1) write(unit=iwrite,FMT=203) iter

write(unit=*, FMT=203) iter

do i=1,Dim_XC

if (method(3)==1) write(unit=iwrite, FMT=202) i, bestmem_XC(i)

write(*,FMT=202) i,bestmem_XC(i)

end do

if (method(3)==1) write(unit=iwrite, FMT=201) bestval

write(unit=*, FMT=201) bestval

end if

iter=iter+1

if ( bestval <= VTR .and. refresh >0) then

write(unit=iwrite, FMT=*) ' The best fitness is smaller than VTR'

write(unit=*, FMT=*) 'The best fitness is smaller than VTR'

exit

endif

end do

!!------end the evolutionary computation------------------------------!!

201 format(2x, 'bestval =', ES14.7, /)

202 format(5x, 'bestmem_XC(', I3, ') =', ES12.5)

203 format(2x, 'No. of iteration =', I8)

end subroutine DE_Fortran90

function randperm(num)

use data_type, only : IB, RP

implicit none

integer(kind=IB), intent(in) :: num

integer(kind=IB) :: number, i, j, k

integer(kind=IB), dimension(num) :: randperm

real(kind=RP), dimension(num) :: rand2

intrinsic random_number

call random_number(rand2)

do i=1,num

number=1

do j=1,num

if (rand2(i) >rand2(j)) then

number=number+1

end if

end do

do k=1,i-1

if (rand2(i) <= rand2(k) .and. rand2(i) >= rand2(k)) then

number=number+1

end if

end do

randperm(i)=number

end do

return


欢迎分享,转载请注明来源:内存溢出

原文地址: http://outofmemory.cn/yw/12546825.html

(0)
打赏 微信扫一扫 微信扫一扫 支付宝扫一扫 支付宝扫一扫
上一篇 2023-05-26
下一篇 2023-05-26

发表评论

登录后才能评论

评论列表(0条)

保存