IF statement formatting (related to Runge-Kutta ODE solver)

Help
1

Hello All,
I have inherited this block of fortran code however, it doesn’t compile. I don’t like using code I don’t understand, but did want to figure this one out.
It seems like the “if” statement is formatted in some abbreviated form. Is this legitimate? To me it appears to be a “If-Then-Else” statement. Anyone know what this is? It is the fourth line down.

PROGRAM RUNGECALL
DIMENSION :: Y(10),F(10)
READ (*,*) X,N,XLIM,H,M,(Y(I),I=1,N)
8 IF(X-XLIM)6,6,7
6 CALL RUNGE(N,Y,F,X,H,M,K)
GO TO (10,20),K
10 F(1)=Y(2)*Y(2)-1-Y(1)*Y(3)
F(2)=Y(1)
F(3)= Y(2)
GO TO 6
20 WRITE (*,*) X,(Y(I),i=1,N)
GO TO 8
7 STOP
END PROGRAM RUNGECALL
2

Line 4 is an arithmetic if, and what it means and how it should be replaced is described for example here.

I wonder why the code does not compile for you. When I save it with the .f90 extension indicating free source form and compile, I get

c:\fortran\test>gfortran -c xrungecall.f90
xrungecall.f90:4:17:

    4 | 8 IF(X-XLIM)6,6,7
      |                 1
Warning: Fortran 2018 deleted feature: Arithmetic IF statement at (1)

So there is an appropriate warning, but it can be compiled by gfortran and other compilers.

Running into arithmetic if is exhibit 101 in what puts newcomers off FORTRAN …

3 Likes
3

This problem was solved numerically by Hiemenz in 1911. This was just a couple of years following the development of boundary layer theory by Ludwig Prandtl (1904) and the first application of the theory to the simple case of flow over a flat plate by Paul Richard Heinrich Blasius (1907). Amazingly, a report of Hiemenz’s work is available. (It is, as may be expected, in German. Google Translate does a good job with the text, but smashes all the equations.)

Here is a modern Fortran version of your code (after correction) :

      program rungecall
      implicit none
      real f(10), h, x, xlim, y(10)
      integer i, k, m, n
!
      read(*,*)x, n, xlim, h, m, (y(i),i=1,n)
      do while ( x <= xlim )
         do
            call runge(n,y,f,x,h,m,k)
            if ( k == 2 ) exit
            f(1) = y(2)*y(2) - 1 - y(1)*y(3)
            f(2) = y(1)
            f(3) = y(2)
         end do
         write(*,*)x, (y(i),i=1,n)
      end do
      end program rungecall
4 Likes
4

It appears that the runge subroutine uses a reverse communication interface, quoting Victor Eijkhout:

The idea is to make the function data-independent by keeping
the data entirely outside its scope. Thus for every operation on that data
you have to return, with a request parameter indicating what
operation has to be performed outside, and you have to have a way
to indicate where you resume after the operation is done.

I’ve written code like this, and it’s kinda neat, but I’m not
entirely unconvinced that it’s not totally evil. Or not.

With your subroutine, depending upon the value of k the calling routine must either perform output or evaluate the right hand side function f. By the looks of things x (independent variable) is probably incremented within the runge routine and h is the (variable?) step size. The purpose of m remains unclear, perhaps it’s the interval at which to perform output?

1 Like
5

Hi Guys,
Wow, much Thanks - you guys are great.

Ivanpribec,
m is an index used in the subroutine which must be set to zero before the first call.

Beliavsky,
exhibit 101 what puts newcomers off, LOL. Yes it’s particularly fustrating - thank you for clearing that up.
I do get the same warning, but there were other errors in the code preventing it from running.
Rick

PS I apparently don’t have the Runge subroutine to call, so it looks like I need to write it myself. Just to be clear this is not a pre canned subroutine I can download from somewhere is it?

6

It’s kind of odd to find a program using go to, line numbers, and arithmetic if formatted as a free-form file. I suspect your file initially looked like this:

Click here

rungecall.for:

      PROGRAM RUNGECALL
      DIMENSION :: Y(10),F(10)
      READ (*,*) X,N,XLIM,H,M,(Y(I),I=1,N)
    8 IF(X-XLIM)6,6,7
    6 CALL RUNGE(N,Y,F,X,H,M,K)
      GO TO (10,20),K
   10 F(1)=Y(2)*Y(2)-1-Y(1)*Y(3)
      F(2)=Y(1)
      F(3)= Y(2)
      GO TO 6
   20 WRITE (*,*) X,(Y(I),i=1,N)
      GO TO 8
    7 STOP
      END PROGRAM RUNGECALL

I haven’t found any routine called runge with a Netlib search. But it shouldn’t be too hard to replace it with a different ODE solver.

If you are not bound to Fortran for performance or other reasons, and you’d just like to know the solution to the ODE system:

\begin{bmatrix} \dot{y}_1 \\[1ex] \dot{y}_2 \\[1ex] \dot{y}_3 \end{bmatrix} = \begin{bmatrix} y_2^2 - 1 - y_1 y_3 \\[1ex] y_1 \\[1ex] y_2 \end{bmatrix}

it might be easier to just use Scipy or MATLAB. If I’m not mistaken this system is equivalent to the third order equation f''' + f'' f = (f')^2 - 1.

1 Like
7

Btw, my goal is to make Fortran as easy as MATLAB or SciPy, so that one can do this within Fortran easily.

8

Regardless of which language or software is used, there is more to this problem than meets the eye, and the OP will benefit from doing some reading on this old problem (Falkner Skan equation from fluid dynamics).

The Runge-Kutta method is primarily an algorithm for IVPs (Initial Value Problem). The present problem, however, is a BVP (boundary value problem).

One method for solving this BVP is to solve a sequence of IVPs with different values of f’’(0), aiming to satisfy the condition f’(infinity) = 1. This is the so called “shooting method”.

There are other algorithms for BVPs, as well.

[Added on 8 November 2021] I used the BVP solver of Shampine, Muir and Xu and found the solution f’’(0) = 1.23259, which agrees with the solution given in page 156, Table 3.4 of https://engineering.purdue.edu/~aerodyn/AAE416/Spring%202010/StagPt/NotesFromWhite.pdf .

2 Likes
9

A number of ODE solvers are in ode , and more in Collected Algorithms of the ACM .

2 Likes
10

Hello All,
Lots of good information here and I cannot tell you how thankful I am. Thank you. A lot of it is a little over my head, but you guys are really good. You correctly deduced what I am working on from the limited info I have given you.

I am a graduate student studying viscous fluid flows. My undergraduate was back in 1991 so it has been some time and it was Fortran 77 when I worked on projects last.

So Ivanpribec,
I am tasked with performing the numerical analysis preferably using Fortran.

Mecej4,
You are correct this problem is Hiemenz flow - stagnation point, and it is complicated because we do not know y3(0) and I am not clear if I need to include that functionality into our program. The way the problem is worded we are allowed to “guess” at the value and see how well it satisfies y2=1 at infinity.
He stated our code should not be more than 100 lines max so I guess some advanced knowledge of the approximate value and just showing a successive convergence on y2=1 should be sufficient to proceed.

2 Likes
11

Hello Ivanpribec,
I did manage to find the old code to the Runge Subroutine. I took a stab at converting it myself. there was a “computed GO TO” which I converted to a Select case construct and the DO loops had to modified. Not sure I did that right, but I have to study the math now to see if the program makes sense. If you can spot anything I did wrong intuitively here I would appreciate it.

Old Code:

Subroutine Runge (n,y,f,x,h,m,k)
!
!This routine performs runge-kutta calculations by gills method.
!
Dimension Y(10), F(10), Q(10)
m=m+1
go to (1,4,5,3,7),m
1 do 2 i=1,n
2 q(i)=0
A=.5
go to 9
!
3 a=1.707107
! if you need more accuracy use a=1.7071067811865475244
!
4 x=x+.5*h
5 do 6 i=1,n
y(i)=y(i)+a*(f(i)*h-Q(i))
6 q(i)=2.*a*h*f(i)+(1.-3.*a)*q(i)
a=.2928932
! if you need more accuracy use a=.2928932188134524756
!
go to 9
7 do 8 i=1,n
8 Y(i)=y(i)+h*F(i)/6.-q(i)/3.
m=0
k=2
go to 10
9 k=1
10 return
end

New code:

! Subroutine Runge(n,y,f,x,h,m,k)

! n = Number of differential equations to be solved (set by the programmer.
! y = Array of n dependent variables (with initial values set by the programmer.
! f = Array of the n derivatives of the variables y (the programmer must provide
! an expression in the main program for the calculation of each F(i))
! x = Independent variable (Initialized by the programmer).
! h = step size delta x (set by the programmer)
! m = Index used in the subroutine which must be set equal to zero by the programmer
! before the first call.
! k = Integer from the subroutine which is used as the argument of a
! computed GO TO Statement (Converted to a do while or Select case Construct here)
!
Subroutine Runge (n,y,f,x,h,m,k)
!
!This subroutine performs runge-kutta calculations by gills method.
!
implicit none
real Y(10), F(10), Q(10),x,h,a
integer n,m,k,i
!
m=m+1
!
Select Case (m)
case (1)
do i=1,n
q(i)=0
A=.5
k=1
End do
!
Case (2)
x=x+.5*h
do i=1,n
y(i)=y(i)+a*(f(i)*h-Q(i))
q(i)=2.*a*h*f(i)+(1.-3.*a)*q(i)
a=.2928932
k=1
! if you need more accuracy use a=.2928932188134524756
end do
!
Case (3)
do i=1,n
y(i)=y(i)+a*(f(i)*h-Q(i))
q(i)=2.*a*h*f(i)+(1.-3.*a)*q(i)
a=.2928932
k=1
! if you need more accuracy use a=.2928932188134524756
end do
!
Case(4)
a=1.707107
! if you need more accuracy use a=1.7071067811865475244
x=x+.5*h
do i=1,n
y(i)=y(i)+a*(f(i)*h-Q(i))
q(i)=2.*a*h*f(i)+(1.-3.*a)*q(i)
a=.2928932
k=1
! if you need more accuracy use a=.2928932188134524756
end do
!
Case (5)
do i=1,n
Y(i)=y(i)+h*F(i)/6.-q(i)/3.
m=0
k=2
end do
End select
!
return
end
12 
Subroutine Runge (n,y,f,x,h,m,k)
!
!This routine performs runge-kutta calculations by gills method.
!
Dimension Y(10), F(10), Q(10)
m=m+1
go to (1,4,5,3,7),m
1 do 2 i=1,n
2 q(i)=0
A=.5
go to 9
!
3 a=1.707107
! if you need more accuracy use a=1.7071067811865475244
!
4 x=x+.5*h
5 do 6 i=1,n
y(i)=y(i)+a*(f(i)*h-Q(i))
6 q(i)=2.*a*h*f(i)+(1.-3.*a)*q(i)
a=.2928932
! if you need more accuracy use a=.2928932188134524756
!
go to 9
7 do 8 i=1,n
8 Y(i)=y(i)+h*F(i)/6.-q(i)/3.
m=0
k=2
go to 10
9 k=1
10 return
end
13

Here is an adaptation of the RK4 code at RosettaCode:

program rungekutta
    implicit none
    integer, parameter :: dp = kind(1d0)
    integer, parameter :: n  = 3 !order of ODE
    real(dp) :: t, dt, tstart, tstop
    real(dp), dimension(n) :: y, k1, k2, k3, k4
 
    tstart = 0.0d0
    tstop =  5.0d0    ! sufficiently large approximation of infinity
    dt = 0.01d0
    y(1:2)  = 0.0d0   ! Initial conditions at t = tstart
    print'(A,$)','Enter guess for f"(0): '
    read *,y(3)       ! Initial condition at t = tstart
    t = tstart
    write (6, '(a,f4.1,a,3f12.8)') 'y(', t, ') = ', y
!    
    do while (t < tstop)
       k1 = dt*f(t, y)
       k2 = dt*f(t+dt/2, y+k1/2)
       k3 = dt*f(t+dt/2, y+k2/2)
       k4 = dt*f(t+dt, y+k3)
       y = y+(k1+2*(k2+k3)+k4)/6
       t = t+dt
       if (abs(nint(t)-t) <= 1d-12) then
          write (6, '(a,f4.1,a,3f12.8)') 'y(', t, ') = ', y
       end if
    end do
!
contains
    function f(t,y)
       real(dp), intent(in) :: t, y(n)
       real(dp) :: f(n)
       f = [y(2),y(3),y(2)*y(2)-1d0-y(1)*y(3)] !Hiemenz equation
    end function f
end program rungekutta

The program prompts for you to feed it a guess for f’’(0). Try different values in the neighborhood of 1.2, “shooting” for the output solution to satisfy f’(5) = 1, f’’(5) = 0. After doing this, you can experiment with altering the upper limit of t to values such as 4 and 6.

1 Like
14

Hello mecej4,
Thank you so much for this. Between all of the responses here with this on I have basically what I need now to proceed on this assignment.
I think I am just going to start a new thread here - as this one has progressed far beyond the original question.
Much thanks guys!
Rick