Okay, I beat my head upon a wall for a while on this one. Why does this work:
y = y+(k1+2*(k2+k3)+k4)/6
But this does not?
y+(1/6)*(k1+(2*k2)+(2*k3)+k4)
Because you have to write it like this:
y+(1._dp/6)*(k1+(2*k2)+(2*k3)+k4)
Where dp is the kind of your (double precision?) number. If it is just single precision, you can write it as:
y+(1./6)*(k1+(2*k2)+(2*k3)+k4)
If you write just 1/6, the answer is 0.
2 Likes
It is important to use informative compiler options, especially when debugging. For the program
integer :: k
k = 3
print*,(1/6)*k
end
gfortran -Wall -Wextra xdiv.f90 says
xdiv.f90:3:9:
print*,(1/6)*k
1
Warning: Integer division truncated to constant '0' at (1) [-Winteger-division]
1 Like
Got it, thank you guys…
RP
Hi Guys,
Here is my first working pass on the program. I think I should add some code to prevent an infinite loop if you input bad values for the guess, but if you put them in right (1.3 and 1.1 for instance). The program works.
If you can be so kind as to point out ways to tighten up the code I’d really appreciate it.
Program AE6009
implicit none
!
integer, parameter :: dp = kind(1d0)
integer, parameter :: n = 3 !Order of ODE's
real(dp) :: a, b, c, g !Variables for shooting method
real(dp) :: eta, deta, etastart, etastop !Variables for Runge-cutta
real(dp), dimension(n) :: y, k1, k2, k3, k4
!
Write (6,*)''
write (6,*) 'Enter a guess for f"(0) < 0'
Read (*,*) a
Write (6,*) 'Enter a guess for f"(0) > 0'
Read (*,*) b
!
g = 1.0d0
do while (g > 1d-10 .or. g < -1d-10)
!
C=0.5*(a+b)
etastart = 0.0d0
etastop = 4.8d0 ! eta at infinity
deta = 0.01d0
y(1) = 0.0d0 ! Initial Condition f(0) = 0
y(2) = 0.0d0 ! Initial Condition f'(0) = 0
y(3) = c ! Initial Condition guess for f''(0)
eta = etastart
!
do while(eta < etastop)
call Runge ()
if (abs(nint(10d0*eta)-(10d0*eta)) <= 1d-10) then
end if
end do
!
g = y(2)-1d0
!
if (g > 0d0)then
a = c
else
b = c
end if
end do
write (6,*)''
Write (6,*) 'The value for f"(0) that satisfies fprime(inf) = 1'
Write (6,'(a,f10.5)') 'f"(0) =' ,c
Write (6,*)''
!
! Now calculate the solutions to the ODE's
!
etastart = 0.0d0
etastop = 4.8d0 ! eta at infinity
deta = 0.01d0
y(1) = 0.0d0 ! Initial Condition f(0) = 0
y(2) = 0.0d0 ! Initial Condition f'(0) = 0
y(3) = c ! Initial Condition f"(0) = c
eta = etastart
Write (6, '(a,f4.1,a,3f10.5)') 'y(',eta,') = ',y
!
do while(eta < etastop)
call Runge ()
if (abs(nint(10d0*eta)-(10d0*eta)) <= 1d-12) then
Write (6, '(a,f4.1,a,3f10.5)') 'y(',eta,') = ',y
end if
end do
!
! Program subroutines and sub functions
!
Contains
Subroutine Runge ()
k1 = deta*f(eta,y) !Runge-cutta 4th order formulas
k2 = deta*f(eta+deta/2, y+k1/2)
k3 = deta*f(eta+deta/2, y+k2/2)
k4 = deta*f(eta+deta, y+k3)
y = y+((1d0/6)*(k1+(2*k2)+(2*k3)+k4))
eta = eta+deta
end subroutine
!
Function f(eta,y)
real(dp), intent(in) :: eta, y(n)
real (dp) :: f(n)
f = [y(2), y(3), y(2)*y(2)-1d0-y(1)*y(3)] ! Hiemenz Equation
end Function
!
end Program AE6009
It looks good, but I have inserted some suggestions and questions as comments prefaced by !!.
Program AE6009
implicit none
!
integer, parameter :: dp = kind(1d0)
integer, parameter :: n = 3 !Order of ODE's
real(dp) :: a, b, c, g !Variables for shooting method
real(dp) :: eta, deta, etastart, etastop !Variables for Runge-cutta
real(dp), dimension(n) :: y, k1, k2, k3, k4
!
Write (6,*)'' !! use * instead of 6 or output_unit from iso_fortran_env
write (6,*) 'Enter a guess for f"(0) < 0'
Read (*,*) a
Write (6,*) 'Enter a guess for f"(0) > 0'
Read (*,*) b
!
g = 1.0d0
do while (g > 1d-10 .or. g < -1d-10)
!
C=0.5*(a+b) !! why is the do while block not indented? Also, use 0.5_dp instead of 0.5
etastart = 0.0d0 !! this group of statements is repeated below. It could be made
!! an internal subroutine
etastop = 4.8d0 ! eta at infinity
deta = 0.01d0
y(1) = 0.0d0 ! Initial Condition f(0) = 0
y(2) = 0.0d0 ! Initial Condition f'(0) = 0
y(3) = c ! Initial Condition guess for f''(0)
eta = etastart
!
do while(eta < etastop)
call Runge ()
if (abs(nint(10d0*eta)-(10d0*eta)) <= 1d-10) then
end if !! should there be a statement before this?
end do
!
g = y(2)-1d0
!
if (g > 0d0)then !! put a space before then for readability
a = c
else !! should have same indentation as if
b = c
end if
end do
write (6,*)''
Write (6,*) 'The value for f"(0) that satisfies fprime(inf) = 1'
Write (6,'(a,f10.5)') 'f"(0) =' ,c
Write (6,*)''
!
! Now calculate the solutions to the ODE's
!
etastart = 0.0d0
etastop = 4.8d0 ! eta at infinity
deta = 0.01d0
y(1) = 0.0d0 ! Initial Condition f(0) = 0
y(2) = 0.0d0 ! Initial Condition f'(0) = 0
y(3) = c ! Initial Condition f"(0) = c
eta = etastart
Write (6, '(a,f4.1,a,3f10.5)') 'y(',eta,') = ',y
!
do while(eta < etastop)
call Runge () !! should have same indentation as line below
if (abs(nint(10d0*eta)-(10d0*eta)) <= 1d-12) then
Write (6, '(a,f4.1,a,3f10.5)') 'y(',eta,') = ',y
end if
end do
!
! Program subroutines and sub functions
!
Contains
Subroutine Runge ()
k1 = deta*f(eta,y) !Runge-cutta 4th order formulas
k2 = deta*f(eta+deta/2, y+k1/2)
k3 = deta*f(eta+deta/2, y+k2/2)
k4 = deta*f(eta+deta, y+k3)
y = y+((1d0/6)*(k1+(2*k2)+(2*k3)+k4))
eta = eta+deta
end subroutine
!
Function f(eta,y)
real(dp), intent(in) :: eta, y(n)
real (dp) :: f(n)
f = [y(2), y(3), y(2)*y(2)-1d0-y(1)*y(3)] ! Hiemenz Equation
end Function
!
end Program AE6009
1 Like
Beliavsky,
Thank you very much for your comments. I’ll work on them all.
Is there any way I can reduce or consolidate this output into perhaps one line?
Write (output_unit,*)''
Write (output_unit,*) "The value for f''(0) that satisfies f'(inf) = 1 is:"
Write (output_unit,'(a,f10.5)') ' f"(0) =' ,c
Write (output_unit,*)''
Yes, use the slash / edit descriptor, described for example here. For example, first two lines above can be written
Write (output_unit,"(/,a)") "The value for f''(0) that satisfies f'(inf) = 1 is:"
Fortran formats have lots of obscure but useful features. I wonder what the best guide to learning them is.
Hi Guys,
Okay, some typo’s aside, I think it’s done. Any other recommendations to tighten it up? I did notice that if I enter guesses too high I get NaN errors in the run. I am not sure how to avoid that.
RP
Program AE6009
use iso_fortran_env !To specify common output_unit
implicit none !To declare all vars explicitly
!
integer, parameter :: dp = kind(1d0)
integer, parameter :: n = 3 !Order of ODE's
real(dp) :: a, b, c, g !Variables for shooting method
real(dp) :: eta, deta, etastart, etastop !Variables for Runge-cutta
real(dp), dimension(n) :: y, k1, k2, k3, k4
integer :: i
!
Write (output_unit,*)''
Write (output_unit,*) "Enter a guess for f''(0) > 0"
Read (*,*) a
Write (output_unit,*) "Enter a guess for f''(0) < 0"
Read (*,*) b
!
! Test for valid input variables
!
If (a < 0 .or. b<0) then
Write (output_unit,*) 'Negative values for guesses not permitted'
Stop
End if
!
CAll Initialize_Variables()
y(3) = a
Do while(eta < etastop)
Call Runge ()
End do
If (Y(2) < 0d0) then
Write (output_unit,*) "Bad guess for f''(0) > 0"
Write (output_unit,'(a,f10.5)') ' f"(0) = ' ,y(2)
Stop
End if
!
Call Initialize_Variables()
Y(3) = b
Do while(eta < etastop)
Call Runge ()
End do
If (Y(2) > 0d0) then
Write (output_unit,*) "Bad guess for f''(0) < 0"
Write (output_unit,'(a,f10.5)') ' f"(0) = ' ,y(2)
Stop
End if
!
g = 1.0d0 !We need some IC for g
i = 0
Do while (g > 1d-10 .or. g < -1d-10)
i = i+1
C = 0.5d0 * (a + b)
Call Initialize_Variables()
Do while(eta < etastop)
Call Runge ()
If (abs(nint(10d0*eta)-(10d0*eta)) <= 1d-10) then
End if
End do
g = y(2)-1d0
If (g > 0d0)then
a = c
Else
b = c
End if
If (i > 1d4) then
Write (output_unit,*) "Input guesses too far apart, infinite loop aborted"
Stop
End if
End do
Write (output_unit,"(/a)") "The value for f''(0) that satisfies f'(inf) = 1 is:"
Write (output_unit,'(a,f10.5/a)') ' f"(0) =' ,c
!
! Now calculate the solutions to the ODE's
!
Call Initialize_Variables()
Write (output_unit, '(a,f4.1,a,3f10.5)') 'y(',eta,') = ',y
!
Do while(eta < etastop)
Call Runge ()
If (abs(nint(10d0*eta)-(10d0*eta)) <= 1d-12) then
Write (output_unit, '(a,f4.1,a,3f10.5)') 'y(',eta,') = ',y
End if
End do
!
! Program subroutines and sub functions
!
Contains
Subroutine Initialize_Variables()
etastart = 0.0d0
etastop = 4.8d0 ! eta at infinity
deta = 0.01d0
y(1) = 0.0d0 ! Initial Condition f(0) = 0
y(2) = 0.0d0 ! Initial Condition f'(0) = 0
y(3) = c ! Initial Condition f"(0) = c
eta = etastart
End subroutine
!
Subroutine Runge ()
k1 = deta*f(eta,y) !Runge-cutta 4th order formulas
k2 = deta*f(eta+deta/2, y+k1/2)
k3 = deta*f(eta+deta/2, y+k2/2)
k4 = deta*f(eta+deta, y+k3)
y = y+((1d0/6)*(k1+(2*k2)+(2*k3)+k4))
eta = eta+deta
End subroutine
!
Function f(eta,y)
real(dp), intent(in) :: eta, y(n)
real (dp) :: f(n)
f = [y(2), y(3), y(2)*y(2)-1d0-y(1)*y(3)] ! Hiemenz Equation
End Function
!
End Program AE6009
Id like to output the data to a CSV file or something I can import into a charting software. Is there a simple way to do this?
Rick
Something similar to
write (unit, 100) eta, y
100 format(F4.1,3(',',F10.5))
should suffice, assuming that you take care of opening the file before writing to it, and closing the file after all the writing has been done.
Hello All,
Just to update you, I finished the program and submitted it. It works great! With all of your help it really made it much easier for me. Thank you.
Only thing I didn’t resolve which was not necessary, if entered a value for f"(0) which was too large I got NaN errors in the output file. I assume those are negative infinity errors due to DP numbers getting so large, but not really sure. I just instruct the user to input closer values.
If you have any comments on the NaN errors I would appreciate it.
If there are any questions about the program please let me know. FORTRAN is a very powerful tool, I’m glad I was forced to revisit it. Looking forward to another project I can use it on…
Kind regards,
Rick
Glad it worked for you. Note that since the 1990 standard, the official spelling is Fortran, with only the first letter capitalized. Wikipedia says
The names of earlier versions of the language through FORTRAN 77 were conventionally spelled in all-uppercase (FORTRAN 77 was the last version in which the Fortran character set included only uppercase letters. The official language standards for Fortran have referred to the language as “Fortran” with (rather than “FORTRAN” in all-uppercase) since Fortran 90.
1 Like
Hi Beliavsky,
Thanks, as if I didn’t feel old enough already. 
The “kids” in my graduate class weren’t even born yet when I graduated with my UG. It’s funny to hear them all complain about having to program in a machine type language like this. None of them get it. I couldn’t even imagine them having to run a program from a terminal prompt instead of icon on the desk top.
It was FORTRAN 77 I studied in 1989.
Thanks for the help.
Rick
@RickP330 , are the other students in your class heavily into non-STEM pursuits for some reason? Or at commercial companies influenced by Windows OS or something?
Because otherwise with the increasing popularity of interpreted and interactive computing environments (Python, R. etc.) with the so-called REPL mode and also Linux (so much work from command line), what people are doing now is conceptually not all that different from what I recently saw in a video at the enterprise where I work of late 1970s of employees doing calculations using various job scripts on terminals connected to an enterprise mainframe system.