Differential equation for radioactive cooling in fortran

Question

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

  END SUBROUTINE dydx

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
    
  END SUBROUTINE heun

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)
ENDDO

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:

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.

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;
    
            Heun(t_0,x_0,h,N);
    
        }// 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]);
            }

            try
            {
                // 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
                arquivo0.close();
            }
            catch (Exception e)
            {
                System.out.println("Erro "+e.getMessage());
            }
            finally
            {
                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.