Intel® Fortran Compiler
Build applications that can scale for the future with optimized code designed for Intel® Xeon® and compatible processors.

Hull Code in Fortran 66

JohnNichols
Valued Contributor III
759 Views

I stumbled across some old Fortran code - 1966. This is the DVERK routine used to solve ODE's, in this case first order. 

It was developed by Hull out of Uni of Toronto in 76. 

I stumbled across a reference in looking at the old Conte and de Boore book on Numerical Analysis - I used that in 1978 in a numerical analysis class. I have a copy on my desktop - it is useful for theory. 

I realize this is an MKL type problem - but what the heck. 

My interest is in looking to solve the stability problems in the 2nd order RKN solver -  . 

John

0 Kudos
6 Replies
JohnNichols
Valued Contributor III
759 Views
program ODE

    implicit none
    Integer N, NW, IND, K
    double precision X, XEND,Y(2), TOL,C(24),W(2,9)
    external fcn1

    N=2
    X =0.0
    IND = 1
    NW = 2
    TOL = 0.000001
    Y(1) = 1.0
    Y(2) = 3.0

    write(*,*)X,Y(1),Y(2)

    do 10 K = 1,10
        XEND = DFLOAT(K)
        CALL dverk(N,FCN1,X,Y,XEND,TOL,IND,C,NW,W)
        write(*,*)X,Y(1),Y(2)
10  end do


    end program ODE

   subroutine FCN1(N,X,Y,YP) 
       Implicit None
       Integer N
       Double Precision X, Y(N), YP(N)
      YP(1) = 2.00*Y(1) *(1.0 - Y(2))
      YP(2) = Y(2)*(Y(1) - 1.0)
      return
      end subroutine FCN1

Why do i need the external FCN1 in the main program -- I do not routinely pass subroutines as arguments 

DVERK running is attached -- only required a few implicits and the calling routines 

answers differ in the 10th place from the publlished results 

John 

0 Kudos
mecej4
Honored Contributor III
759 Views

Nichols, John wrote:
Why do i need the external FCN1 in the main program -- I do not routinely pass subroutines as arguments 

Procedures (functions or subroutines) need to be passed as arguments when solving a problem where the problem description contains an arbitrary nonlinear function. Examples:

  • Solve F(x) = 0
  • Integrate F(x) from x = a to x = b
  • Find the minimum point of F(x)
  • Fit y(x) to a set of data pairs (x1, y1), ... , (xn,yn)

More complex examples involve F being a function of many independent variables.

If you use IMPLICIT NONE or imply it with a compiler option, and provide no declaration for FCN1, it would be flagged as missing type information. Otherwise, it would be a local variable of type REAL, with undefined value; this is not what we wanted!

From another point of view, EXTERNAL is needed to distinguish data variables from code variables.

If we provide explicit interfaces, EXTERNAL is not needed.

0 Kudos
mecej4
Honored Contributor III
759 Views

Verner has an interesting Web page on R-K-V methods at http://people.math.sfu.ca/~jverner/ .

0 Kudos
JohnNichols
Valued Contributor III
759 Views

Thanks for the explanation. The program wokrs nicely, the details of the algorithms at the start of the program are great and the code is like a spaghetti western -- My name is noboby springs to mind.

Runge Kutta are fascinating. The problem is the rather interesting problem of adding the velocity component to the second order problem.  I have the RKN method for the general equation sort of working, but the issue is uncontrolled growth on the second derivative -- I am now looking through these other versions to work out how to control it -- my method is to simple at this stage.  

I was playing with the RKN code yesterday for a small bridge and the results showed once the truck passed over the bridge the next couple of seconds showed higher displacement in a sine wave, Jim D commented on this point elsewhere and he was correct.

 I was somewhat taken aback by the result wondering what I had done wrong when I observed in my recent development of the accel velocity and displacement plots for the bridge from experimental data the same thing. 

The truck is such a heavy load relative to the bridge that the dynamics is impacted by the loss of the truck load. 

 

0 Kudos
FortranFan
Honored Contributor II
759 Views

Code appears closer to FORTRAN 77 than 66; intrinsic functions used in this code such as DMAX1 (the use of the generic MAX would be what the current standard suggests) were not available in FORTRAN 66.

0 Kudos
JohnNichols
Valued Contributor III
759 Views

I think it has been played with over the years -- it is still quite jumbled.   I was going over the code in pencil the other night -- take a bit to unwind it. 

0 Kudos
Reply