program main
c*********************************************************************72
c
cc STRING_SIMULATION() solves and plots a solution of the 1d wave equation.
c
c Discussion:
c
c This program solves the 1D wave equation of the form:
c
c Utt = c^2 Uxx
c
c over the spatial interval (X1,X2) and time interval (T1,T2),
c with initial conditions:
c
c U(T1,X) = U_T1(X),
c Ut(T1,X) = UT_T1(X),
c
c and boundary conditions of Dirichlet type:
c
c U(T,X1) = U_X1(T),
c U(T,X2) = U_X2(T).
c
c The value C represents the propagation speed of waves.
c
c The program uses the finite difference method, and marches
c forward in time, solving for all the values of U at the next
c time step by using the values known at the previous two time steps.
c
c Central differences may be used to approximate both the time
c and space derivatives in the original differential equation.
c
c Thus, assuming we have available the approximated values of U
c at the current and previous times, we may write a discretized
c version of the wave equation as follows:
c
c Uxx(T,X) = ( U(T, X+dX) - 2 U(T,X) + U(T, X-dX) ) / dX^2
c Utt(T,X) = ( U(T+dt,X ) - 2 U(T,X) + U(T-dt,X ) ) / dT^2
c
c If we multiply the first term by C^2 and solve for the single
c unknown value U(T+dt,X), we have:
c
c U(T+dT,X) = ( C^2 * dT^2 / dX^2 ) * U(T, X+dX)
c + 2 * ( 1 - C^2 * dT^2 / dX^2 ) * U(T, X )
c + ( C^2 * dT^2 / dX^2 ) * U(T, X-dX)
c - U(T-dT,X )
c
c (Equation to advance from time T to time T+dT, except for FIRST step)
c
c However, on the very first step, we only have the values of U
c for the initial time, but not for the previous time step.
c In that case, we use the initial condition information for dUdT
c which can be approximated by a central difference that involves
c U(T+dT,X) and U(T-dT,X):
c
c dU/dT(T,X) = ( U(T+dT,X) - U(T-dT,X) ) / ( 2 * dT )
c
c and so we can estimate U(T-dT,X) as
c
c U(T-dT,X) = U(T+dT,X) - 2 * dT * dU/dT(T,X)
c
c If we replace the "missing" value of U(T-dT,X) by the known values
c on the right hand side, we now have U(T+dT,X) on both sides of the
c equation, so we have to rearrange to get the formula we use
c for just the first time step:
c
c U(T+dT,X) = 1/2 * ( C^2 * dT^2 / dX^2 ) * U(T, X+dX)
c + ( 1 - C^2 * dT^2 / dX^2 ) * U(T, X )
c + 1/2 * ( C^2 * dT^2 / dX^2 ) * U(T, X-dX)
c + dT * dU/dT(T, X )
c
c (Equation to advance from time T to time T+dT for FIRST step.)
c
c It should be clear now that the quantity ALPHA = C * DT / DX will affect
c the stability of the calculation. If it is greater than 1, then
c the middle coefficient 1-C^2 DT^2 / DX^2 is negative, and the
c sum of the magnitudes of the three coefficients becomes unbounded.
c
c Licensing:
c
c This code is distributed under the MIT license.
c
c Modified:
c
c 25 December 2012
c
c Author:
c
c John Burkardt
c
c Local Parameters:
c
c Local, double precision ALPHA, the CFL stability parameter.
c
c Local, double precision C, the wave speed.
c
c Local, double precision DT, the time step.
c
c Local, double precision DX, the spatial step.
c
c Local, integer M, the number of time steps.
c
c Local, integer N, the number of spatial intervals.
c
c Local, double precision T1, T2, the initial and final times.
c
c Local, double precision U(M+1,N+1), the computed solution.
c
c Local, double precision X1, X2, the left and right spatial endpoints.
c
implicit none
integer m
parameter ( m = 30 )
integer n
parameter ( n = 40 )
double precision alpha
double precision c
integer command_unit
integer data_unit
double precision dt
double precision dx
double precision f
double precision g
integer i
integer j
double precision k
double precision t
double precision t1
double precision t2
double precision u(0:m,0:n)
double precision x
double precision x1
double precision x2
c = 0.25D+00
t1 = 0.0D+00
t2 = 3.0D+00
x1 = 0.0D+00
x2 = 1.0D+00
call timestamp ( )
write ( *, '(a)' ) ''
write ( *, '(a)' ) 'string_simulation():'
write ( *, '(a)' ) ' Fortran77 version'
write ( *, '(a)' )
& ' Simulate the behavior of a vibrating string.'
dx = ( x2 - x1 ) / real ( n, kind = 8 )
dt = ( t2 - t1 ) / real ( m, kind = 8 )
alpha = ( c * dt / dx ) ** 2
write ( *, '(a,g14.6)' ) ' ALPHA = ( C * dT / dX )^2 = ', alpha
c
c Warn the user if ALPHA will cause an unstable computation.
c
if ( 1.0D+00 .lt. alpha ) then
write ( *, '(a)' ) ''
write ( *, '(a)' ) ' Warning!'
write ( *, '(a)' ) ' ALPHA is greater than 1.'
write ( *, '(a)' ) ' The computation is unstable.'
end if
c
c Time step 0:
c Use the initial condition for U.
c
u(0,0) = 0.0D+00
do j = 1, n - 1
x = dble ( j ) * dx
u(0,j) = f ( x )
end do
u(0,n) = 0.0D+00
c
c Time step 1:
c Use the initial condition for dUdT.
c
u(1,0) = 0.0D+00
do j = 1, n - 1
x = dble ( j ) * dx
u(1,j) =
& ( alpha / 2.0D+00 ) * u(0,j-1)
& + ( 1.0D+00 - alpha ) * u(0,j)
& + ( alpha / 2.0D+00 ) * u(0,j+1)
& + dt * g ( x )
end do
u(1,n) = 0.0D+00
c
c Time steps 2 through M:
c
do i = 2, m
u(i,0) = 0.0D+00
do j = 1, n - 1
u(i,j) =
& alpha * u(i-1,j-1)
& + 2.0D+00 * ( 1.0D+00 - alpha ) * u(i-1,j)
& + alpha * u(i-1,j+1)
& - u(i-2,j)
end do
u(i,n) = 0.0D+00
end do
c
c Write data file.
c
call get_unit ( data_unit )
open ( unit = data_unit, file = 'string_data.txt',
& status = 'replace' )
do i = 0, m
t = dble ( i ) * dt
do j = 0, n
x = dble ( j ) * dx
write ( data_unit, '(f10.4,2x,f10.4,2x,f10.4)' )
& x, t, u(i,j)
end do
write ( data_unit, '(a)' ) ''
end do
close ( unit = data_unit )
write ( *, '(a)' ) ''
write ( *, '(a)' )
& ' Plot data written to the file "string_data.txt".'
c
c Write gnuplot command file.
c
call get_unit ( command_unit )
open ( unit = command_unit, file = 'string_commands.txt' )
write ( command_unit, '(a)' ) 'set term png'
write ( command_unit, '(a)' ) 'set output "string.png"'
write ( command_unit, '(a)' ) 'set grid'
write ( command_unit, '(a)' ) 'set style data lines'
write ( command_unit, '(a)' ) 'unset key'
write ( command_unit, '(a)' ) 'set xlabel "<---X--->"'
write ( command_unit, '(a)' ) 'set ylabel "<---Time--->"'
write ( command_unit, '(a,i4,a)' )
& 'splot "string_data.txt" using 1:2:3 with lines'
write ( command_unit, '(a)' ) 'quit'
close ( unit = command_unit )
write ( *, '(a)' )
& ' Gnuplot commands written to the file "string_commands.txt".'
c
c Terminate.
c
write ( *, '(a)' ) ''
write ( *, '(a)' ) 'STRING_SIMULATION:'
write ( *, '(a)' ) ' Normal end of execution.'
write ( *, '(a)' ) ''
call timestamp ( )
stop
end
function f ( x )
c*********************************************************************72
c
cc F supplies the initial condition.
c
c Licensing:
c
c This code is distributed under the MIT license.
c
c Modified:
c
c 01 December 2012
c
c Author:
c
c John Burkardt
c
c Parameters:
c
c Input, double precision X, the location.
c
c Output, double precision F, the value of the solution at time 0 and location X.
c
implicit none
double precision f
double precision x
if ( 0.25D+00 <= x .and. x <= 0.50D+00 ) then
f = ( x - 0.25D+00 ) * ( 0.50D+00 - x )
else
f = 0.0D+00
end if
return
end
function g ( x )
c*********************************************************************72
c
cc G supplies the initial derivative.
c
c Licensing:
c
c This code is distributed under the MIT license.
c
c Modified:
c
c 01 December 2012
c
c Author:
c
c John Burkardt
c
c Parameters:
c
c Input, double precision X, the location.
c
c Output, double precision G, the value of the time derivative of the solution
c at time 0 and location X.
c
implicit none
double precision g
double precision x
g = 0.0D+00
return
end
subroutine get_unit ( iunit )
c*********************************************************************72
c
cc GET_UNIT returns a free Fortran unit number.
c
c Discussion:
c
c A "free" Fortran unit number is a value between 1 and 99 which
c is not currently associated with an I/O device. A free Fortran unit
c number is needed in order to open a file with the OPEN command.
c
c If IUNIT = 0, then no free Fortran unit could be found, although
c all 99 units were checked (except for units 5, 6 and 9, which
c are commonly reserved for console I/O).
c
c Otherwise, IUNIT is a value between 1 and 99, representing a
c free Fortran unit. Note that GET_UNIT assumes that units 5 and 6
c are special, and will never return those values.
c
c Licensing:
c
c This code is distributed under the MIT license.
c
c Modified:
c
c 02 September 2013
c
c Author:
c
c John Burkardt
c
c Parameters:
c
c Output, integer IUNIT, the free unit number.
c
implicit none
integer i
integer iunit
logical value
iunit = 0
do i = 1, 99
if ( i .ne. 5 .and. i .ne. 6 .and. i .ne. 9 ) then
inquire ( unit = i, opened = value, err = 10 )
if ( .not. value ) then
iunit = i
return
end if
end if
10 continue
end do
return
end
subroutine timestamp ( )
c*********************************************************************72
c
cc TIMESTAMP prints out the current YMDHMS date as a timestamp.
c
c Discussion:
c
c This Fortran77 version is made available for cases where the
c Fortran90 version cannot be used.
c
c Licensing:
c
c This code is distributed under the MIT license.
c
c Modified:
c
c 12 January 2007
c
c Author:
c
c John Burkardt
c
c Parameters:
c
c None
c
implicit none
character * ( 8 ) ampm
integer d
character * ( 8 ) date
integer h
integer m
integer mm
character * ( 9 ) month(12)
integer n
integer s
character * ( 10 ) time
integer y
save month
data month /
& 'January ', 'February ', 'March ', 'April ',
& 'May ', 'June ', 'July ', 'August ',
& 'September', 'October ', 'November ', 'December ' /
call date_and_time ( date, time )
read ( date, '(i4,i2,i2)' ) y, m, d
read ( time, '(i2,i2,i2,1x,i3)' ) h, n, s, mm
if ( h .lt. 12 ) then
ampm = 'AM'
else if ( h .eq. 12 ) then
if ( n .eq. 0 .and. s .eq. 0 ) then
ampm = 'Noon'
else
ampm = 'PM'
end if
else
h = h - 12
if ( h .lt. 12 ) then
ampm = 'PM'
else if ( h .eq. 12 ) then
if ( n .eq. 0 .and. s .eq. 0 ) then
ampm = 'Midnight'
else
ampm = 'AM'
end if
end if
end if
write ( *,
& '(i2,1x,a,1x,i4,2x,i2,a1,i2.2,a1,i2.2,a1,i3.3,1x,a)' )
& d, month(m), y, h, ':', n, ':', s, '.', mm, ampm
return
end