✨ 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

Recommended Posts

✨ Easy with Python ☆ Estimated elapsed time after death ✨
[Co-occurrence analysis] Easy co-occurrence analysis with Python! [Python]
Easy folder synchronization with Python
Execution time measurement with Python With
Easy Python compilation with NUITKA-Utilities
Easy HTTP server with Python
Time synchronization (Windows) with Python
Easy time series prediction with Prophet
[Python] Easy parallel processing with Joblib
Easy Python + OpenCV programming with Canopy
Easy email sending with haste python3
Bayesian optimization very easy with Python
Easy data visualization with Python seaborn.
Easy parallel execution with python subprocess
Easy modeling with Blender and Python
[Python] Super easy test with assert statement
[Python] Easy argument type check with dataclass
Just print Python elapsed time in seconds
[Easy Python] Reading Excel files with openpyxl
Easy web app with Python + Flask + Heroku
Easy web scraping with Python and Ruby
[Python] Easy Reinforcement Learning (DQN) with Keras-RL
[Python] Easy introduction to machine learning with python (SVM)
Csv output from Google search with [Python]! 【Easy】
[Python3] A story stuck with time zone conversion
Use logger with Python for the time being
How to measure execution time with Python Part 1
Easy Lasso regression analysis with Python (no theory)
Make your Python environment "easy" with VS Code
How to measure execution time with Python Part 2
Statistics with python
[Super easy] Simultaneous face recognition and facial expression recognition in real time with Python and OpenCV!
Python with Go
Twilio with Python
Integrate with Python
python time measurement
Play with 2016-Python
AES256 with python
Tested with Python
python starts with ()
First time python
with syntax (Python)
Bingo with python
Zundokokiyoshi with python
Microcomputer with Python
Cast with python
Stop EC2 for specified time + start with Lambda (python)
Run with CentOS7 + Apache2.4 + Python3.6 for the time being
Easy partial download of mp4 with python and youtube-dl!
Get standard output in real time with Python subprocess
Easy way to scrape with python using Google Colab
How to measure mp3 file playback time with python
To avoid spending time coloring shapes drawn with python
[Time series with plotly] Dynamic visualization with plotly [python, stock price]