如何使用龙格-库塔法求解牛顿冷却定律 (RK4)

Question

我需要帮助来实施解决牛顿冷却定律的四阶龙格库塔方法。由于这个问题引入了时间 (t)，我对给定条件下的位置感到困惑。这是给定的：时间间隔从 t = 0 到 t = 20（以秒为单位），对象温度 = 300，环境温度 = 70，时间增量为 .1，比例常数 = 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);

}

}

下面是我求解 ODE 的 RK4 方法。在下面的代码中，假设 y(1) = 10 和 delta x = 0.001

，我将 ODE y' = y - x 求解为近似 y(1.005)

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

}

}

根据公式 T(t) = Ts + (T0 - Ts) * e^(-k*t) 我应该有 75.1 的近似值来求解牛顿方程。 Ts = 环境温度，T0 = 物体初始温度，t = 20（经过的秒数），k = .19 比例常数

Answer 1

我猜（但不是很难）您要求解的 ODE 是

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

如您所见，右侧并不直接取决于时间。由于您没有尝试为系统编写代码，因此环境温度 T_amb 很可能是一个常数。因此移动常量并使用一致的函数名称并将 ODE 函数参数返回格式 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);

    }


}

如何使用龙格-库塔法求解牛顿冷却定律 (RK4)

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

java

math

for-loop

runge-kutta