HDQ Basith Studio

Specialized for graphic design and 3D modelling.

Download anything

Here you can get many useful software, music, and many kind of file.

Get windows problem?

Don't worry, here you may get the solutions of your windows problems.

Tips and Tricks

Sometimes a problem can be solved with a simple solution. Find the easiest and strangest solutions here.

Wanna get a high degree?

يرفع الله الذين آمنوا منكم والذين أوتوا العلم درجات

Minggu, 10 Juni 2012

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

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

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

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

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

Senin, 04 Juni 2012

Buletin Embun edisi 2