Today I'm trying to evaluate this differential equation for internal energy in a gas in Fortran:

$$ \frac{du}{dt} = - \frac{n_H^2}{\rho}\frac{\Lambda}{n_H^2} $$

Where nH is the density of hydrogen in the gas (costant), u is the internal energy: $$ u(T) = \frac{1}{\gamma - 1}\frac{kT}{\mu m_p} $$

And lambda is the cooling function (depending also only on temperature), sum of all the radioactive processes in the gas that make it cool down over time.

I'm trying to solve this using a RK2 method in Fortran (Heun), and I'm sure that the actual solving algorithm works fine, but I'm not sure on how to even implement this.

I'm using this subroutine to define the right hand side of the equation:

  SUBROUTINE dydx(neq, y, f)
    INTEGER, INTENT(IN) :: neq
    REAL*8, DIMENSION(neq), INTENT(IN) :: y
    REAL*8, DIMENSION(neq), INTENT(OUT) :: f
    f= -y


Already, I'm pretty sure this is not the right way to do this. I'm trying to make it in a way such that $\frac{dy}{dx} = -f(x,y)$, but I really don't get how.

This is my Heun alg:

SUBROUTINE heun(neq, h, yold, ynew)
    INTEGER, INTENT(IN) :: neq
    REAL*8, INTENT(IN) :: h 
    REAL*8, DIMENSION(neq), INTENT(IN) ::yold
    REAL*8, DIMENSION(neq), INTENT(OUT) :: ynew

    REAL*8, DIMENSION(neq) :: f, ftilde
    INTEGER :: i

    CALL dydxv(neq, yold, f)

    DO i=1, neq
      ynew(i) = yold(i) + h*f(i)
    END DO

    CALL dydxv(neq, ynew, ftilde)

    DO i=1, neq
      ynew(i) = yold(i) + 0.5d0*h*(f(i) + ftilde(i))
    END DO

And I think this one's fine. This is how I'm calling it (lambda, E_0-the constants-, density are arrays of size n where I stored all the values given by their respective functions in the temperatures I'm considering):

DO i = 1,n
  CALL heun(n, h, -lambda/(E_0*density), temperature)

I want to plot the results against a timescale in units of cooling time, which I can calculate analytically as: $$ t_{cool} = \frac{u}{(\frac{\Lambda}{\rho})} $$

And this is the graph I get using this implementation:

enter image description here

Which looks off but at least it's going down, since it's supposed to be cooling down.

My final goal is to have a starting temperature of a million Kelvin and simulate cooling until ten thousand Kelvin, and see how long the process takes (compared to the analytical cooling time). How do I implement this boundary conditions?

I know I'm asking many questions but I hardly know where to start when it comes to numerical differentiation. If you could at least give me a starting point or good resources on where to get started, it would be very much appreciated.

  • 2
    $\begingroup$ Is the first equation even right? The $n_H^2$ terms cancel, for example. $\endgroup$ Commented Sep 5, 2022 at 22:18
  • 3
    $\begingroup$ It's also unclear what the role of $u(t)$ is. In the first equation, it appears to be the unknown you want to solve for. But then in the next equation, it appears that you already know what $u$ is (and that it depends on some variable $T$ rather than $t$). Can you clarify? $\endgroup$ Commented Sep 6, 2022 at 19:55
  • $\begingroup$ Yes they cancel, I kept them there to remind myself that I have to calculate the cooking function per unit density. I want to solve for T (temperature) but the problem is that I don’t have an explicit function of T(u), which would make it easier to solve that. I now figured out that I can find this function using a secant method on the function of T that I get from u(T). Lowercase t is time, which neither function directly depends on $\endgroup$ Commented Sep 6, 2022 at 20:44
  • $\begingroup$ Could you provide the various constants that are needed shown in the formulation, along with the cooling function lambda. I will then code a solution in Python. $\endgroup$
    – Tarik
    Commented Nov 28, 2023 at 12:26
  • $\begingroup$ Could you also please unlearn Real*8 which is not Fortran and has never been part of Fortran, may not be supported by your compiler, and may not do what you think it does. See stackoverflow.com/questions/838310/fortran-90-kind-parameter $\endgroup$
    – Ian Bush
    Commented Nov 28, 2023 at 16:17

1 Answer 1


Well, I didn't understand your specific physical problem very well. But about the numerical method of Heun, I made a code in Java so that you can take it as a base if you want to rewrite it in Fortran. In this code, we generate a txt (or dat) file with the results, and we can use the gnuplot to generate the graph.

    import javax.swing.*;
    import java.io.FileWriter;
     * Author: Carlos Eduardo da Silva Lima
     * Modified Euler method (Heun)
     * Enter t_0: start time
     * Enter x_0: starting position
     * Enter h: Step h, from Heun's algorithm
     * Enter N: Iteration amount
     * Textbook/Biography: Numerical Calculus Neide Berthold Franco
    public class Heun
        public static void main(String[] args)
            double t_0 = 0.0;
            double x_0 = 1.0;
            double h   = 1E-4;
            int    N   = 100000;
        }// End of main method main
        // Enter the expression for the derivative dy(t,x)/dt = f(t,x)
        public static double f(double t, double x)
            return -1.2*x+7*Math.exp(-0.3*t);
        // This method implements Heun's algorithm
        public static void Heun(double t_0, double x_0, double h, int N)
            double k1, k2, k3;
            double [] t = new double[N]; // time t
            double [] x = new double[N]; // position in the inactive time t, x(t)

            t[0] = t_0; // Initial time t0
            x[0] = x_0; // Initial position x0, in the initial time t0
            for(int i = 0; i<=N-2; i++)
                k1 = f(t[i],x[i]);
                k2 = f(t[i]+(h/3),x[i]+((h*k1)/3));
                k3 = f(t[i]+((2*h)/3),x[i]+((2*h*k2)/3));
                x[i+1] = x[i] + ((h*(k1+3*k3))/4);
                t[i+1] = t[i] + h;
            // On-screen output (Promp)
            for(int i = 0; i<=N-2; i++)
                // System.out.printf("t = %f | x = %f\n",t[i],x[i]);

                // Creation of .txt or .dat file
                FileWriter arquivo0 = new FileWriter("Heun.txt",false);
                for(int j = 0; j<(N-1); j++)
                    arquivo0.write(t[j]+"       "+x[j]+"       "+"\n");
                // Closing the .txt or .dat file
            catch (Exception e)
                System.out.println("Erro "+e.getMessage());
                System.out.println("Fim do bloco try-catch-finally");
                JOptionPane.showMessageDialog(null,".txt(or .dat) files created successfully!");
        }// End of Heun's method
    }// End of main class Heun

I also add the excerpt from the "Cálculo Numérico" textbook, page 407 with the equations for the Heun method.

enter image description here


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.