Search code examples
javafor-loopmathrunge-kutta

How to solve Newton's Cooling Law using Runge-Kutta method (RK4)


I need help implementing my working 4th-Order Runge-Kutta method for solving Newton's Cooling Law. Since time (t) is introduced with this problem, I am confused on the placings for the given conditions. Here is what's given: time interval begins t = 0 to t = 20 (in seconds), object temp = 300, ambient temp = 70, time increment is .1, and constant of proportionality = 0.19

public class RungeKutta {

public static double NewtonsCoolingLaw(double objectTemp,double  ambientTemp)
     {
         double k = 0.19;       
         return -k * (objectTemp - ambientTemp);
     }

public static void main(String[] args) {

    double result = 0.0;  
    double  initialObjectTemp = 300.0, givenAmbientTemp = 70.0;
    double deltaX = (20.0 - 0)/10000;

    for(double t = 0.0; t <= 20.0; t += .1)
    {        
        double k1 = deltaX * NewtonsLaw(initialObjectTemp,givenAmbientTemp);
        double k2 = deltaX * NewtonsLaw(initialObjectTemp + (deltaX/2.0),givenAmbientTemp + (k1/2.0));
        double k3 = deltaX * NewtonsLaw(initialObjectTemp + (deltaX/2.0), givenAmbientTemp + (k2/2.0));
        double k4 = deltaX * NewtonsLaw(initialObjectTemp + deltaX, givenAmbientTemp + k3);
        givenAmbientTemp = givenAmbientTemp + (1.0/6.0) * (k1 + (2.0 * k2) + (2.0 * k3) + k4);
        result = givenAmbientTemp;
    }       
    System.out.println("The approx. object temp after 20 seconds is: " + result);

}

}

Bellow is my RK4 method for solving ODEs. In the code below, I solve the ODE y' = y - x to approximate y(1.005) given that y(1) = 10 and delta x = 0.001

public class RungeKutta {

public static double functionXnYn(double x,double  y)
     {
         return y-x;
     }

public static void main(String[] args) {

    double deltaX = (1.005 - 0)/10000;
    double y = 10.0;
    double result = 0.0;  

    for(double x = 1.0; x <= 1.005; x = x + deltaX)
    {
        double k1 = deltaX * functionXnYn(x,y);
        double k2 = deltaX * functionXnYn(x + (deltaX/2.0),y + (k1/2.0));
        double k3 = deltaX * functionXnYn(x + (deltaX/2.0), y + (k2/2.0));
        double k4 = deltaX * functionXnYn(x + deltaX, y + k3);
        y = y + (1.0/6.0) * (k1 + (2.0 * k2) + (2.0 * k3) + k4);
        result = y;
    }       
    System.out.println("The value of y(1.005) is: " + result);

}

}

Based on the formula T(t) = Ts + (T0 - Ts) * e^(-k*t) I should have an approximation of 75.1 for solving Newton's DE. Ts = ambient temp, T0 = object initial temp, t = 20 (seconds elapsed), and k = .19 constant of proportionality


Solution

  • I'm guessing (but not really hard) that the ODE you are trying to solve is

    dT(t)/dt = -k*(T(t)-T_amb)
    

    As you can see, the right side does not directly depend on the time. As you make no attempts to code for a system, it is likely that the ambient temperature T_amb is a constant. Thus moving the constants around and using consistent function names and returning the ODE function arguments to the format time, state variable

    public class RungeKutta {
    
        public static double CoolingLaw(double time, double objectTemp)
             {
                 double k = 0.19, ambientTemp = 70.0;     
                 return -k * (objectTemp - ambientTemp);
             }
    
        public static void main(String[] args) {
    
            double result = 0.0;  
            double  objectTemp = 300.0;
            double dt = 0.1
    
            for(double t = 0.0; t <= 20.0; t += dt)
            {        
                double k1 = dt * CoolingLaw(t, objectTemp);
                double k2 = dt * CoolingLaw(t + (dt/2.0), objectTemp + (k1/2.0));
                double k3 = dt * CoolingLaw(t + (dt/2.0), objectTemp + (k2/2.0));
                double k4 = dt * CoolingLaw(t + dt, objectTemp + k3);
                objectTemp = objectTemp + (1.0/6.0) * (k1 + (2.0 * k2) + (2.0 * k3) + k4);
                result = objectTemp;
            }       
            System.out.println("The approx. object temp after 20 seconds is: " + result);
    
        }
    
    
    }