Contoh Runge-Kutta dengan Fortran

 program main
!====================================================================
! Initial value problem for a single 1st-order ODE
! Call Methods: Euler, Predictor-Corrector, Runge-Kutta 4th order
!====================================================================
implicit none
double precision dx, dt, xi, ti, xf, tf, tmax
integer i, key
external dx

! open files (if needed)
!open (unit=6, file='result1.dat')

! initial conditions
xi = 1.00
ti = 0.00
! "time" step and max time
dt =   0.1
tmax = 2.0
! select a method:
! key=1 Euler, key=2 predictor-corrector, key=3 RK 4th order
key = 3

!* print the header and initial conditions
write (*,*) ' 1st order single ODE '
if (key==1) write (*,*) ' Simple Euler method  '
if (key==2) write (*,*) ' Predictor-Corrector  '
if (key==3) write (*,*) ' Runge-Kutta 4th order'
write (*,*) '     t         x(t)   '
write (6,100) ti, xi


do while(ti <= tmax)
  tf = ti + dt

if (key==1)  call euler1(dx,ti,xi,tf,xf)
if (key==2)  call euler1m(dx,ti,xi,tf,xf)
if (key==3)  call RK4d11(dx,ti,xi,tf,xf)

  write (6,100) tf, xf
  ti = tf
  xi = xf
end do

100 format(2f10.5)
end program main

  function dx(t,x)
!----------------------------------------------
! dx(t,x)- function dx/dt in the 1st order ODE
!----------------------------------------------
implicit none
double precision dx, x, t
  dx = (-1.0)*x     ! solution x=exp(-t)
! the carrying capacity problem
!  dx = 0.1*x-0.0002*x**2
end function dx

  subroutine RK4d11(dx,ti,xi,tf,xf)
!==================================================
! Solution of a single 1st order ODE dx/dt=f(x,t)
! Method: 4th-order Runge-Kutta method
! Alex G. February 2010
!--------------------------------------------------
! input ...
! dx(t,x)- function dx/dt (supplied by a user)
! ti - initial time
! xi  - initial position
! tf  - time for a solution
! output ...
! xf  - solution at point tf
!==================================================
implicit none
double precision dx,ti,xi,tf,xf
double precision h,k1,k2,k3,k4

h  = tf-ti

k1 = h*dx(ti,xi)
k2 = h*dx(ti+h/2.0,xi+k1/2.0)
k3 = h*dx(ti+h/2.0,xi+k2/2.0)
k4 = h*dx(ti+h,xi+k3)

xf = xi + (k1 + 2.0*(k2+k3) + k4)/6.0

end subroutine RK4d11

  subroutine euler1m(dx,ti,xi,tf,xf)
!==================================================
! Solution of a single 1st order ODE dx/dt=f(x,t)
! Method:  Modified Euler (predictor-corrector)
! Alex G. February 2010
!--------------------------------------------------
! input ...
! dx(t,x)- function dx/dt (supplied by a user)
! ti - initial time
! xi  - initial position
! tf  - time for a solution
! output ...
! xf  - solution at point tf
!==================================================
implicit none
double precision dx,xi,ti,xf,tf

xf = xi + dx(ti,xi)*(tf-ti)
xf = xi + (dx(ti,xi)+dx(tf,xf))*0.5*(tf-ti)

end subroutine euler1m

  subroutine euler1(dx,ti,xi,tf,xf)
!==================================================
! Solution of a single 1st order ODE dx/dt=f(x,t)
! Method:  Simple Euler (predictor-corrector)
! Alex G. February 2010
!--------------------------------------------------
! input ...
! dx(t,x)- function dx/dt (supplied by a user)
! ti - initial time
! xi  - initial position
! tf  - time for a solution
! output ...
! xf  - solution at point tf
!==================================================
implicit none
double precision dx,xi,ti,xf,tf

xf = xi + dx(ti,xi)*(tf-ti)

end subroutine euler1

Selengkapnya...

Contoh Differensial dengan Fortran

!====================================================================
!  Driver program: Calculates first- or second-order derivatives
!  for a function defined in ni base points
!  Method: based on explicit differentiation of Lagrange interpolation
!  AG: October 2009
!====================================================================
implicit none
integer, parameter :: ni=101    ! initial arrays size for interpolation
integer, parameter :: nc=21  ! number of points where derivatives to be calc.
double precision xmin, xmax  ! interval
double precision xi(ni), yi(ni)
double precision x, yx, yxx, step, f, deriv3
double precision fx, fxx
integer i, j


! open files
!open (unit=1, file='table01.dat')
!open (unit=2, file='table02.dat')

xmin =  0.0
xmax =  3.1415926

!
!   step 1: generate xi and yi from f(x) xmin, xmax, nint
!
step = (xmax-xmin)/float(ni-1)
do i=1,ni
  xi(i) = xmin + step*float(i-1)
  yi(i) = f(xi(i))
!  write (*,200) xi(i), yi(i)
end do

!
!  step 2: calculate derivatives
!
step = (xmax-xmin)/float(nc-1)
write (*,202)
do i=1, nc
  x = xmin + step*float(i-1)
  yx = cos(x)
  fx = deriv3(x, xi, yi, ni, 1)
  yxx = -sin(x)
  fxx= deriv3(x, xi, yi, ni, 2)
  write (*,201) x, yx, fx, yxx, fxx
end do

200 format (2f12.5)
201 format (7f12.5)
202 format ('       x           fx         fx_num', &
            '      fxx        fxx_num')
!203 format (' Orders of divided difference interpolation')

stop
end

!
!  Function f(x)
!
  function f(x)
  implicit none
  double precision f, x
  f = sin(x)
  return
  end

  function deriv3(xx, xi, yi, ni, m)
!====================================================================
! Evaluate first- or second-order derivatives
! using three-point Lagrange interpolation
! written by: Alex Godunov (October 2009)
!--------------------------------------------------------------------
! input ...
! xx    - the abscissa at which the interpolation is to be evaluated
! xi()  - the arrays of data abscissas
! yi()  - the arrays of data ordinates
! ni - size of the arrays xi() and yi()
! m  - order of a derivative (1 or 2)
! output ...
! deriv3  - interpolated value
!============================================================================*/

implicit none
integer, parameter :: n=3
double precision deriv3, xx
integer ni, m
double precision xi(ni), yi(ni)
double precision x(n), f(n)
integer i, j, k, ix

! exit if too high-order derivative was needed,
if (m > 2) then
  deriv3 = 0.0
  return
end if

! if x is ouside the xi(1)-xi(ni) interval set deriv3=0.0
if (xx < xi(1) .or. xx > xi(ni)) then
  deriv3 = 0.0
  return
end if

! a binary (bisectional) search to find i so that xi(i-1) < x < xi(i)
i = 1
j = ni
do while (j > i+1)
  k = (i+j)/2
  if (xx < xi(k)) then
    j = k
  else
    i = k
  end if
end do

! shift i that will correspond to n-th order of interpolation
! the search point will be in the middle in x_i, x_i+1, x_i+2 ...
  i = i + 1 - n/2

! check boundaries: if i is ouside of the range [1, ... n] -> shift i
if (i < 1) i=1
if (i + n > ni) i=ni-n+1

!  old output to test i
!  write(*,100) xx, i
!  100 format (f10.5, I5)

! just wanted to use index i
ix = i

! initialization of f(n) and x(n)
do i=1,n
  f(i) = yi(ix+i-1)
  x(i) = xi(ix+i-1)
end do

! calculate the first-order derivative using Lagrange interpolation
if (m == 1) then
    deriv3 =          (2.0*xx - (x(2)+x(3)))*f(1)/((x(1)-x(2))*(x(1)-x(3)))
    deriv3 = deriv3 + (2.0*xx - (x(1)+x(3)))*f(2)/((x(2)-x(1))*(x(2)-x(3)))
    deriv3 = deriv3 + (2.0*xx - (x(1)+x(2)))*f(3)/((x(3)-x(1))*(x(3)-x(2)))
! calculate the second-order derivative using Lagrange interpolation
  else
    deriv3 =          2.0*f(1)/((x(1)-x(2))*(x(1)-x(3)))
    deriv3 = deriv3 + 2.0*f(2)/((x(2)-x(1))*(x(2)-x(3)))
    deriv3 = deriv3 + 2.0*f(3)/((x(3)-x(1))*(x(3)-x(2)))
end if
end function deriv3

Selengkapnya...

Contoh Integral dengan Fortran

http://www.odu.edu/~agodunov/computing/programs/index.html

    Program main
!=============================================
! Integration of a function using Simpson rule
!=============================================
implicit none
double precision f, a, b, integral
integer n, i
double precision, parameter:: pi = 3.1415926
external f

a = 0.0
b = pi

n = 2

write(*,100)

do i=1,16
   call simpson(f,a,b,integral,n)
   write (*,101) n, integral
   n = n*2
end do

100   format('     nint   Simpson')
101   format(i9,1pe15.6)
end


  Function f(x)
!----------------------------------------
! Function for integration
!----------------------------------------
implicit none
double precision f, x
 f = sin(x)
! f = x*cos(10.0*x**2)/(x**2 + 1.0)
return
end

 Subroutine simpson(f,a,b,integral,n)
!==========================================================
! Integration of f(x) on [a,b]
! Method: Simpson rule for n intervals 
! written by: Alex Godunov (October 2009)
!----------------------------------------------------------
! IN:
! f   - Function to integrate (supplied by a user)
! a      - Lower limit of integration
! b      - Upper limit of integration
! n   - number of intervals
! OUT:
! integral - Result of integration
!==========================================================
implicit none
double precision f, a, b, integral,s
double precision h, x
integer nint
integer n, i

! if n is odd we add +1 to make it even
if((n/2)*2.ne.n) n=n+1

! loop over n (number of intervals)
s = 0.0
h = (b-a)/dfloat(n)
do i=2, n-2, 2
   x   = a+dfloat(i)*h
   s = s + 2.0*f(x) + 4.0*f(x+h)
end do
integral = (s + f(a) + f(b) + 4.0*f(a+h))*h/3.0
return
end subroutine simpson

Selengkapnya...

Contoh Metode Eliminasi Gauss dengan Fortran

PROGRAM GAUSSIAN

PARAMETER (IN=20)

REAL:: A(IN,IN), B(IN)

PRINT *,'INPUT NUMBER OF EQUATIONS N='

READ (5,*) N

PRINT *,'INPUT MATRIX COEFFICIENTS A(I,J)='

READ (5,*) ((A(I,J),J=1,N),I=1,N)

PRINT *,'INPUT RIGHT-HAND SIDE VECTOR B(I)='

READ (5,*) (B(I),I=1,N)

WRITE (6,*) ('****GAUSSIAN ELIMINATION ****')

WRITE (6,*)

WRITE (6,*) ('COEFFICIENT MATRIX you inputed:')

CALL PRINTA(A,IN,N,N,6)

WRITE(6,*)

WRITE(6,*) ('right hand side vector you inputed:')

CALL PRINTV(B,N,6)

WRITE(6,*)

! CONVERT TO UPPER TRIANGULAR FORM

DO K = 1,N-1

IF (ABS(A(K,K)).GT.1.E-6) THEN

DO I = K+1, N

X = A(I,K)/A(K,K)

DO J = K+1, N

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

ENDDO

B(I) = B(I) - B(K)*X

ENDDO

ELSE

WRITE (6,*) 'ZERO PIVOT FOUND IN LINE:'

WRITE (6,*) K

STOP

END IF

ENDDO

WRITE(6,*) 'MODIFIED MATRIX'

CALL PRINTA(A,IN,N,N,6)

WRITE(6,*)

WRITE(6,*) 'MODIFIED RIGHT HAND SIDE VECTOR'

CALL PRINTV (B,N,6)

WRITE(6,*)

! BACK SUBSTITUTION

DO I = N,1,-1

SUM = B(I)

IF (I.LT.N) THEN

DO J= I+1,N

SUM = SUM - A(I,J)*B(J)

ENDDO

END IF

B(I) = SUM/A(I,I)

ENDDO

! PRINT THE RESULTS

write(6,*) ('SOLUTION VECTOR')

CALL PRINTV(B,N,6)

!

END PROGRAM GAUSSIAN

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



SUBROUTINE PRINTA(A,IA,M,N,ICH)

!

! WRITE A 2D ARRAY TO OUTPUT CHANNEL 'ICH'

!

REAL A(IA,*)

DO I =1,M

WRITE(ICH,2) (A(I,J),J=1,N)

ENDDO

2 FORMAT(1X,6E12.4)

!

END SUBROUTINE PRINTA

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



SUBROUTINE PRINTV(VEC,N,ICH)

!

! WRITE A COLUMN VECTOR TO CHANNEL 'ICH'

!

REAL VEC(*)

WRITE(ICH,1) (VEC(I),I=1,N)

1 FORMAT(1X,6E12.4)

!

END SUBROUTINE PRINTV

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

Selengkapnya...

Buletin Embun edisi 2

Download Buletin Embun Edisi 2

Selengkapnya...