✨ Easy with Python ☆ Estimated elapsed time after death ✨

What if there was a corpse in front of you?

(引用;名探偵コナン「甘く冷たい宅配便」) When Conan finds his body, he thinks about the elapsed time after death at the very beginning, right? If you know the estimated time of death, you can narrow down the number of suspects. \ (^ O ^) / So this time, I will introduce how to find out more about the elapsed time after death!

example

There was a man's body in a refrigerator car set at room temperature of 16 ° C. The man weighed 86 kg, and when Conan measured the rectal temperature of the body, it was 27 ° C. Find the elapsed time after the death of a man.

Newton's law of cooling

First, consider a simplified phenomenon. Consider the problem that hot coffee in a cup is cooled by the ambient temperature (air). At this time, it is said that the following "Newton's law of cooling" holds approximately. The differential equation that expresses Newton's law of cooling is

\frac{dT}{dt} = -γ(T - T_s) (1)

It is represented by. Here, each character is $ T $: Hot water temperature $ T_s $: Ambient temperature $ \ gamma $: Cooling constant Represents. This differential equation can be solved analytically. First, divide both sides of equation (1) by $ T-T_s $. \frac{1}{T-T_s}dT = -γdt

If both sides are integrated and the constant of integration is C, log (T - T_s) = -γt + C It will be. Solving this equation for T gives the following equation. T = T_s + e^{-γt + C} (2)
When the initial condition $ T (0) = T_0 $ is imposed, the constant of integration C becomes C = log (T_0-T_s)

It will be. Substituting this into Eq. (2), the analytical solution is given below. T(t) = T_s + (T_0 - T_s) e^{-γt}

(Also, transform this) \frac{T(t)-T_s}{T_0-T_s} = e^{-γt}

Estimating the elapsed time after death by continuous measurement of rectal temperature

Next, use Newton's law of cooling to determine the elapsed time after death. According to "Estimation of elapsed time after death by continuous measurement of rectal temperature (Waseda University Advanced Science and Engineering)", in the case of a corpse, the differential equation is transformed as follows. \frac{T_r-T_e}{T_o-T_e}=1.25e^{B*t}-0.25e^{5B*t} (T_e\le23.2) \frac{T_r-T_e}{T_o-T_e}=1.11e^{B*t}-0.11e^{10B*t} (T_e\ge23.2) B=-1.2815(bW)^{-0.625}+0.0284 However, $ T_r $ is rectal temperature ℃ $ T_e $ is the ambient temperature ℃ $ T_o $ is the rectal temperature at the time of death and is 37.2 ° C. $ T_e $ is constant at the mean temperature of measurement time at ambient temperature ℃ $ W $ Corpse weight kg $ b = 1 $ is the correction factor $ B $ is Newton's cooling constant $ t $ is the elapsed time after death (h) will do. Let's solve this example immediately based on this!

import sympy

def death_time(Tr,Te,W):
    Tr=Tr
    Te=Te
    W=W
    B = sympy.Symbol('B')
    t = sympy.Symbol('t')
    To = 37.2
    b = 1
    E = sympy.S.Exp1
    equation1 = 1.25*(E**(B*t))-0.25*(E**(5*B*t))-(Tr-Te)/(To-Te)
    equation2 = 1.11*(E**(B*t))-0.11*(E**(10*B*t))-(Tr-Te)/(To-Te)
    equation3 = -1.2815*((b*W)**(-0.625))+0.0284-B
    if Te <= 23.2:
        print(sympy.solve([equation1,equation3]))
    else:
        print(sympy.solve([equation2,equation3]))    
death_time(27,16,86) #Setting this issue

According to the $$ reference, $ T_e $ uses the average ambient temperature, but this time we set the room temperature of 16 ° C at the time of discovery to $ T_e $. You can tell the result without saying it, but scipy couldn't directly calculate including such a power variable, and I continued to calculate even when I woke up in the morning (crying). I decided to solve it based on the idea of trial and error method according to the reference.

(引用KrylovSubspaceAcceleratedNewtonAlgorithm:ApplicationtoDynamicProgressiveCollapseSimulationofFrames)

Since the process of finding the differential of a function and solving the equation each time it is repeated is large, we decided to calculate using the Krilov partial exchange method, which approximates the derivative by the difference as shown in the above figure.

import math
from scipy.optimize import newton_krylov
import sympy

#(3)Corresponds to the formula. body weight(kg)Calculate Newton's cooling constant B from. To(Rectal temperature at death)37.2
def predict_B(W):
    W=W
    B = sympy.Symbol('B')
    b = 1
    equation3 = -1.2815*((b*W)**(-0.625))+0.0284-B
    B = sympy.solve([equation3])[B]
    return B  

#(1)(2)Equivalent to an expression. Case classification by environmental temperature Te
def predict_deathtime(W,Tr,Te,ta):
    Bx = predict_B(W)
    Tr,Te,ta = Tr,Te,ta
    E = sympy.S.Exp1
    To = 37.2
    def F(x):
        if Te <= 23.2:
            return 1.25*(E**(Bx*x))-0.25*(E**(5*Bx*x))-(Tr-Te)/(To-Te)              
        else:
            return 1.11*(E**(Bx*x))-0.11*(E**(10*Bx*x))-(Tr-Te)/(To-Te)
    guess = ta
    sol = newton_krylov(F, guess, method='lgmres')
    print(sol)

#(body weight,Rectal temperature,Environmental temperature,Initial value of elapsed time after death)Enter in the order of
#Note that the differential value cannot be obtained by setting the initial value of the elapsed time after death to 0.
predict_deathtime(86,27,16,2)

Output result 17.19191007875492

Therefore, in this example, it turned out that he died about 17 hours ago!

References

[Let's start computational science! Newton's Law of Cooling ①](https://dekfractal.com/2018/09/23/%E8%A8%88%E7%AE%97%E7%A7%91%E5%AD%A6%E3%82% 92% E3% 81% AF% E3% 81% 98% E3% 82% 81% E3% 82% 88% E3% 81% 86% EF% BC% 81% E3% 80% 80% E3% 83% 8B% E3% 83% A5% E3% 83% BC% E3% 83% 88% E3% 83% B3% E3% 81% AE% E5% 86% B7% E5% 8D% B4% E6% B3% 95% E5% 89% 87 /) Estimation of elapsed time after death by continuous measurement of rectal temperature (Waseda University Advanced Science and Technology) ○ (Study) Mikiyasu Inoue * (Correct) Kiyotaka Sakai (Waseda Univ.) (Correct) Kenichiro Yamamoto (National Defense Medical College) (Correct) Jun Kim