Hull Code in Fortran 66

[JohnNichols]

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

[JohnNichols]

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 

[mecej4]

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 (x₁, y₁), ... , (x_n,y_n)

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.

[mecej4]

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

[JohnNichols]

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. 

[FortranFan]

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.

[JohnNichols]

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.