龙格库塔方法在实际生活中的应用（数值计算Java）

• 提出实际问题

应用于刑事侦查中死亡时间的鉴定，最近一直在看一些经典的悬疑片，里面都会有警察来判断死者的死亡时间，就特别好奇是什么原理，去了解了一下，判断死亡时间是根据牛顿冷却定律推算出来的，那么我就打算用数值分析的方法来实现一下。

• 基本思想或原理

牛顿冷却定律运用了常微分方程的知识，牛顿冷却定律在实际生活中可以解决好多实际问题。牛顿冷却定律的内容是：物体在空气中冷却的速度与物体温度和空气温度之差成正比。

人体体温受大脑中枢神经调节，人死亡后体温调节的功能消失，尸体的温度受外界环境温度的影响。那么我们已知尸体温度变化率服从牛顿冷却定律，即尸体温度的变化率与他周围的温度差成正比。

若物体在房间里，与物体相比，房间非常大，那么我们假设房间的温度为常温，即20摄氏度。物体对房间的温度改变忽略不计。那就可以直接将牛顿冷却定律应用于刑事侦查中死亡时间的鉴定。

当谋杀发生后，尸体温度从原来的36.5摄氏度按照牛顿冷却定律开始下降。周围空气温度始终保持20摄氏度不变，三小时后，尸体温度变为33.5摄氏度。如果尸体被发现时的温度为30摄氏度，假设现在的时间是下午的五点整，那么我们就可以大概知道谋杀是什么时候发生的。由题我们可以得出相应公式：

 由分离变量得：

代入初始条件，可以求得C=16.5，所以得到初值问题的解为:

因为三个小时后尸体的温度为33.5摄氏度，那么我们可以求出

 

• 运行环境

Idea/Java编程语言

• 程序代码

import java.util.ArrayList;
import java.util.List;
/**
 * @author Sun
 * @version 2021.2
 * @date 2022/5/21 10:42
 */
public class Runge_Kutta {
    public static void main(String[] args) {
        List<List<Double>> result = getEstimates(0.0,36.5,0.2,40);
        List<Double> xi = result.get(0);
        List<Double> yi = result.get(1);
        List<Double> reals = result.get(2);
        List<Double> miscount = result.get(3);
        for (int i = 0; i < xi.size(); i++ ) {
  System.out.println("x"+i+"="+xi.get(i)+"\t"+"y"+i+"="+yi.get(i)+"\t"+"real="+reals.get(i)+"\t"+"mis="+miscount.get(i));
        }
    }
    public static List<List<Double>> getEstimates(double x0,double y0,double h,int N) {
        List<List<Double>> result = new ArrayList<>();
        List<Double> xi = new ArrayList<>();
        List<Double> yi = new ArrayList<>();
        List<Double> reals = new ArrayList<>();//精确值
        List<Double> miscount = new ArrayList<>();//误差 real(n)-y(n)
        double x1,y1,yp,yc,real,mis,K1,K2,K3,K4;
        for (int i = 0 ; i < N ; i++ ) {
            x1 = x0 + h;
            K1 = f(x0,y0);
            K2 = f(x0+h/2,y0+h*K1/2);
            K3 = f(x0+h/2,y0+h*K2/2);
            K4 = f(x1,y0+h*K3);
            y1 = y0 + h*(K1 + 2*K2 + 2*K3 + K4)/6;
            real = getReal(x0);
            mis = Math.abs(real - y1);
            xi.add(x1);
            yi.add(y1);
            reals.add(real);
            miscount.add(mis);
            x0 = x1;
            y0 = y1;
        }
        result.add(xi);
        result.add(yi);
        result.add(reals);
        result.add(miscount);
        return result;
    }
    private static double getReal(double x) {
        //精确值 报告中温度函数公式
        return 20+(16.5*Math.pow(2.718281,-0.06689*x));
    }
    private static double f(double x0, double y0) {
        //拟合出的数值解 y‘
        return -1*0.06689*(y0-20);
    }
}

• 运行过程和运行结果

• 

• 实际问题结果

和使用高数里面的常微分方程得到的结果相同，当H=30代入温度函数中，  t=7.4小时，那么我们可以判断谋杀发生的大概时间下午五点前的7.4个小时，也就是早上9时36分。