ODEs in Fortran

As part of my senior thesis work, I’ve been quickly getting up to speed with both programming in Fortran and various numerical methods for solving different types of equations. This, of course, will all come before I learn the relevant topics in my numerical analysis (MTH 451) class, but I don’t mind too much. This thesis is more important than that math course anyway, plus learning it beforehand means I can focus on my other courses at the end of the semester without stressing too much. Sounds good to me!

Anyway, I wrote three programs over the last few hours to solve the ODE u’ = u [ exact solution: u = exp(t) ] just to make sure I was programming the methods right. After all, I may need to use those methods in my final project, so getting them right now is a big step. In all three programs (that differ only in the “solve” function), the step size needs to be changed in two separate places and the function would need to be changed for different problems, both of which would require a re-compile, but for these programs that isn’t too bad. You’d need to do the same for a change to some sort of input file anyway, and I didn’t want to get around that restriction by having user-inputted values. Check them out!

PROGRAM Forward_Euler
    IMPLICIT None
    REAL, PARAMETER :: tstep = 0.2
    REAL, PARAMETER :: tend = 3.0

    PRINT '(A)', "Forward Euler Method"
    CALL solve(tend)

    CONTAINS

    FUNCTION f(u, t)
        REAL :: u, t, f
        f = u
    END FUNCTION f

    FUNCTION advance(u, t) RESULT(unew)
        REAL :: unew, u, t
        unew = u + tstep * f(u, t)
    END FUNCTION advance

    SUBROUTINE solve(tend)
        REAL :: tend, t = 0.0, u = 1.0
        INTEGER :: k = 0
        DO
            IF ( t > tend ) EXIT
            PRINT '(i2, f5.1, f7.3, f7.3, f7.3)', k, t, u, EXP(t), ABS(EXP(t)-u)
            u = advance(u, t)
            t = t + tstep
            k = k + 1
        END DO
    END SUBROUTINE solve
END PROGRAM Forward_Euler
PROGRAM RungeKutta_2
    IMPLICIT None
    REAL, PARAMETER :: tstep = 0.2
    REAL, PARAMETER :: tend = 3.0

    PRINT '(A)', "2nd Order Runge-Kutta Method"
    CALL solve(tend)

    CONTAINS

    FUNCTION f(u, t)
        REAL :: u, t, f
        f = u
    END FUNCTION

    FUNCTION rungekutta2(u, t) RESULT(unew)
        REAL :: u, t, unew, ustep, tnew
        ustep = u + (tstep / 2.0) * f(u, t)
        tnew = t + tstep / 2.0
        unew = u + tstep * f(ustep, tnew)
    END FUNCTION rungekutta2

    SUBROUTINE solve(tend)
        REAL :: tend, t = 0.0, u = 1.0
        REAL, PARAMETER :: tstep = 0.2
        INTEGER :: k = 0.0
        DO
            IF ( t > tend ) EXIT
            PRINT '(i2, f5.1, f7.3, f7.3, f7.3)', k, t, u, EXP(t), ABS(EXP(t)-u)
            u = rungekutta2(u, t)
            t = t + tstep
            k = k + 1
        END DO
    END SUBROUTINE solve
END PROGRAM RungeKutta_2
PROGRAM RungeKutta_4
    IMPLICIT None
    REAL, PARAMETER :: tstep = 0.2
    REAL, PARAMETER :: tend = 3.0

    PRINT '(A)', "4th Order Runge-Kutta Method"
    CALL solve(tend)

    CONTAINS

    FUNCTION f(u, t)
        REAL :: u, t, f
        f = u
    END FUNCTION

    FUNCTION rungekutta4(u, t) RESULT(unew)
        REAL :: u, t, unew, k1, k2, k3, k4, tnew
        k1 = tstep * f(u, t)
        tnew = t + tstep / 2.0
        k2 = tstep * f(u+0.5*k1, tnew)
        k3 = tstep * f(u+0.5*k2, tnew)
        tnew = t + tstep
        k4 = tstep * f(u+k3, t+tstep)
        unew = u + (1.0/6.0) * (k1 + 2*k2 + 2*k3 + k4)
    END FUNCTION

    SUBROUTINE solve(tend)
        REAL :: tend, t = 0.0, u = 1.0
        REAL, PARAMETER :: tstep = 0.2
        INTEGER :: k = 0
        DO
            IF ( t > tend ) EXIT
            PRINT '(i2, f5.1, f7.3, f7.3, f7.3)', k, t, u, EXP(t), ABS(EXP(t)-u)
            u = rungekutta4(u, t)
            t = t + tstep
            k = k + 1
        END DO
    END SUBROUTINE solve
END PROGRAM RungeKutta_4

Again, the only difference is in the “solve” method. Eventually these will all be wrapped up in a module, along with the other algorithms I’ll need for the final project, but for now three separate files is nice. I could put them in a single file and have the program call each one in turn, or only run a user-selected method, but that’s relatively simple compared to actually getting the coding down and to work.

I am very glad that it works.

Advertisements

6 responses to “ODEs in Fortran

  1. this is awesome man

  2. Uday Gaitonde

    I appreciate your neat procedure. I have been a Fortran programmer for about 40 years (for solving algebraic equs, ODEs, PDEs, …) and know the great advantage of doing one thing at a time, from the simple one to the bigger and complex one.

    • Glad I’m getting something right! I was just introduced to Fortran about a month ago (these are all in Fortran 90), and I love the language. May I ask exactly what you need to use it for, past the general ODE, PDE, etc.?

  3. Uday Gaitonde

    I am a Mechanical Engineer and Fortran has been my primary programming language. My work involves solving problems related to thermal engineering – thermodynamics, fluid mechanics, heat transfer, and related aspects of energy conversion. Whether is a simple task of interpolating in a table, or solving transient conduction equations, or solving the NSE (aka CFD!), the language was Fortran. Have programmed in all versions of Fortran: 66, 77, 90. Have not really graduated to 95 and later, as yet.

    • Thanks for the info! I’ll probably need to use it for similar reasons (just related to stars, most likely) over the next few years. For now, though, I’m teaching myself as I go along, which is difficult for me when I’m also worrying about classes and applying to grad school, so I’m sure there’s tons of things I haven’t even encountered that are old hat for you.

  4. Uday Gaitonde

    One thing you need not worry about is Fortran. Just post your queries here, and I will try to help to the extent I can.