微分方程与拟合

刑事侦查中死亡时间的鉴定模型某地发生一起谋杀案，刑侦人员测得尸体温度为30c，此时是下午4点整，假设该人被谋杀前的体温为37c，被杀两个小时后尸体温度为35c，周围空气的温度为20c，试推断谋杀是何时发生的？网上查牛顿冷却定律（若不知道这个定律，就查尸体如何冷却）一、模型假设与变量说明：1、假设尸体的温度按牛顿冷却定律开始下降，即尸体冷却的速度与尸体温度和空气温度之差成正比。（牛顿冷却定律指出：物体在空气中冷却的速度与物体温度和空气温度之差成正比）；2、假设尸体的最初温度为37c,两个小时后尸体温度为35c，且周围空气的温度保持20c不变；3、假设尸体被发现时的温度为30c，时间是下午4点整；4、假设尸体的温度为T（t）（t从谋杀时计）。二、模型的分析与建立：由于尸体的冷却速度dT/dt与尸体温度T和空气温度之差成正比，设比例系数为k（k0为常数），则有dT/dt=-k(T-20)初始条件为T（0）=37c.当函数是递减函数时，正比例系数前是负号。三、模型求解：1、求通解：分离变量得111,20ln20,2020ktcktcktdTkdtTTktceeeTce两端积分得则有T即:2、求特解：将初始条件为T（0）=37c代入通解，得c=17.于是满足该问题的特解为20.0630.0632017,35c352017,0.063,20173010/178.4().ktkttTekekeeth为求出值,根据两小时后尸体温度为这一条件,有求得于是尸体函数为T将T代入上式有,即得用matlab求解dsolve('DT=-0.063*(T-20)',‘T(0)=37')ans=17/exp((63*t)/1000)+20solve('30=20+17*exp(-63*t/1000)','t')ans=-(1000*log(10/17))/63ans=8.4227结果分析于是可以判定谋杀发生在下午4点尸体被发现前的8.4h，即是在上午7点36分发生的。由于温度受很多因素影响，不一定严格成正比，并且与当时死者身体的穿着有关，衣物不同影响温度的变化，所以死亡时间要参考同类案件来确定。农业生产实验模型在研究农业生产的试验中，为分析某地区土豆产量与化肥的关系，得到了每公顷地的氮肥的施肥量与土豆产量的对应关系。根据表格数据给出土豆产量与氮肥施肥量之间的关系。氮肥量（kg）03467101135202259336404471土豆产量(kg)15.1821.3625.7232.2934.0339.4543.1543.4640.8330.75先画散点图x=[03467101135202259336404471];y=[15.1821.3625.7232.2934.0339.4543.1543.4640.8330.75];plot(x,y,'*')05010015020025030035040045050015202530354045拟合x=[03467101135202259336404471];y=[15.1821.3625.7232.2934.0339.4543.1543.4640.8330.75];a=polyfit(x,y,2)a=-0.00030.197114.7416拟合图形x=[03467101135202259336404471];y=[15.1821.3625.7232.2934.0339.4543.1543.4640.8330.75];a=ployfit(x,y,2)x1=0:471;y1=polyval(a,x1);plot(x,y,'*',x1,y1,'r')0501001502002503003504004505001015202530354045x=[03467101135202259336404471];y=[15.1821.3625.7232.2934.0339.4543.1543.4640.8330.75];x1=0:471;y1=-0.0003*x1.^2+0.1971*x1+14.7416;plot(x,y,'*',x1,y1,'r')050100150200250300350400450500101520253035404550再次拟合（以比较）x=[03467101135202259336404471];y=[15.1821.3625.7232.2934.0339.4543.1543.4640.8330.75];b=polyfit(x,y,3)b=-0.0000-0.00020.164515.7098再次画拟合图形x=[03467101135202259336404471];y=[15.1821.3625.7232.2934.0339.4543.1543.4640.8330.75];a=polyfit(x,y,2);b=polyfit(x,y,3);x1=0:471;y1=polyval(a,x1);y2=polyval(b,x1);plot(x,y,'*',x1,y1,'r',x1,y2,'g')0501001502002503003504004505001015202530354045血药浓度模型通过实验测得一次性快速静脉注射300mg药物后的血药浓度数据见表：求血药浓度随时间的变化规律。t/h0.250.511.523468Y(ug/ml)19.2118.1515.3614.1012.899.327.455.243.01先画散点图x=[0.250.511.523468];y=[19.2118.1515.3614.1012.899.327.455.243.01];plot(x,y,'*')0123456782468101214161820根据散点图可知大致呈指数形态（y=a*exp(-bt)）建立M文件（由于是单调递减，所以指数的系数为负的）functiony=fun66(x,t)y=x(1)*exp(-x(2)*t);%输入主程序t=[0.250.511.523468];p=[19.2118.1515.3614.1012.899.327.455.243.01];t0=[100.4292];x=lsqcurvefit('fun66',t0,t,p)p1=fun66(x,t)取t=1,则p=15.36,再令a=10,(15.36=10*exp(-b))于是b=-ln1.536=0.4292x=20.24130.2420p1=Columns1through819.053217.934815.891114.080212.47579.79457.68944.7394Column92.9211建立M文件:functiony=fun66(x,t)y=x(1)*exp(-x(2)*t);%输入主程序t=[0.250.511.523468];p=[19.2118.1515.3614.1012.899.327.455.243.01];t0=[100.4292];x=nlinfit(t,p,'fun66',t0)p1=fun66(x,t)拟合t=[0.250.511.523468];p=[19.2118.1515.3614.1012.899.327.455.243.01];t1=0.25:0.01:8;p1=20.2413*exp(-0.2420*t1);plot(t,p,'*')holdonplot(t1,p1,'g:')0123456782468101214161820012345678246810121416182022t1=0.25:0.01:8;p1=20.2413*exp(-0.2420*t1);plot(t1,p1,'g:')t1=0.25:0.01:8;p1=fun66(x,t1);plot(t1,p1,'g:')转化为线性拟合来完成)(,ln,ln,lnln)(,2121线性的则该方程变为令两边同时取对数t：Y）baYy（btayeaybtt=[0.250.511.523468];y=[19.2118.1515.3614.1012.899.327.455.243.01];Y=log(y);p=polyfit(t,Y,1)p=-0.23472.9943即te：ybeea2347.09943.2219714.192347.0,9714.192347.0,9943.21所以于是symstt=0.25:0.01:8;y=19.9714*exp(-0.2347*t);plot(t,y,'g:')t0=[0.250.511.523468];p=[19.2118.1515.3614.1012.899.327.455.243.01];holdonplot(t0,p,'*')0123456782468101214161820饮酒驾车的拟合图t=[0.250.50.7511.522.533.544.55678910111213141516];p=[3068758282776868585150413835282518151210774];t1=0.25:0.01:16;p1=114.4326*(exp(-0.1855*t1)-exp(-2.0079*t1));plot(t,p,'*')holdonplot(t1,p1,'g:')02468101214160102030405060708090