微分方程与拟合

1、刑事侦查中死亡时间的鉴定模型,某地发生一起谋杀案，刑侦人员测得尸体温度为30c，此时是下午4点整，假设该人被谋杀前的体温为37c，被杀两个小时后尸体温度为35c，周围空气的温度为20c，试推断谋杀是何时发生的？,网上查牛顿冷却定律（若不知道这个定律，就查尸体如何冷却）,一、模型假设与变量说明： 1、假设尸体的温度按牛顿冷却定律开始下降，即尸体冷却的速度与尸体温度和空气温度之差成正比。（牛顿冷却定律指出：物体在空气中冷却的速度与物体温度和空气温度之差成正比）； 2、假设尸体的最初温度为37c,两个小时后尸体温度为35c，且周围空气的温度保持20c不变； 3、假设尸体被发现时的温度为30c，时间是

2、下午4点整； 4、假设尸体的温度为T（t）（t从谋杀时计）。,二、模型的分析与建立：,由于尸体的冷却速度dT/dt与尸体温度T和空气温度之差成正比，设比例系数为k（k0为常数），则有dT/dt=-k(T-20) 初始条件为T（0）=37c. 当函数是递减函数时，正比例系数前是负号。,三、模型求解：,1、求通解： 分离变量得,2、求特解：将初始条件为T（0）=37c代入通解，得c=17.于是满足该问题的特解为,用matlab求解,dsolve(DT=-0.063*(T-20),T(0)=37) ans = 17/exp(63*t)/1000) + 20 solve(30=20+17*exp(-6

3、3*t/1000),t) ans = -(1000*log(10/17)/63 ans = 8.4227,结果分析,于是可以判定谋杀发生在下午4点尸体被发现前的8.4h，即是在上午7点36分发生的。 由于温度受很多因素影响，不一定严格成正比，并且与当时死者身体的穿着有关，衣物不同影响温度的变化，所以死亡时间要参考同类案件来确定。,农业生产实验模型,在研究农业生产的试验中，为分析某地区土豆产量与化肥的关系，得到了每公顷地的氮肥的施肥量与土豆产量的对应关系。根据表格数据给出土豆产量与氮肥施肥量之间的关系。,先画散点图,x=0 34 67 101 135 202 259 336 404 471 ;

4、y=15.18 21.36 25.72 32.29 34.03 39.45 43.15 43.46 40.83 30.75; plot(x,y,*),拟合,x=0 34 67 101 135 202 259 336 404 471 ; y=15.18 21.36 25.72 32.29 34.03 39.45 43.15 43.46 40.83 30.75; a=polyfit(x,y,2) a = -0.0003 0.1971 14.7416,拟合图形,x=0 34 67 101 135 202 259 336 404 471 ; y=15.18 21.36 25.72 32.29 34.0

5、3 39.45 43.15 43.46 40.83 30.75; a=ployfit(x,y,2) x1=0:471; y1=polyval(a,x1); plot(x,y,*,x1,y1, r),x=0 34 67 101 135 202 259 336 404 471 ; y=15.18 21.36 25.72 32.29 34.03 39.45 43.15 43.46 40.83 30.75; x1=0:471; y1=-0.0003*x1.2+0.1971*x1+14.7416; plot(x,y,*,x1,y1, r),再次拟合（以比较）,x=0 34 67 101 135 202

6、259 336 404 471 ; y=15.18 21.36 25.72 32.29 34.03 39.45 43.15 43.46 40.83 30.75; b=polyfit(x,y,3) b= -0.0000 -0.0002 0.1645 15.7098,再次画拟合图形,x=0 34 67 101 135 202 259 336 404 471 ; y=15.18 21.36 25.72 32.29 34.03 39.45 43.15 43.46 40.83 30.75; a=polyfit(x,y,2); b=polyfit(x,y,3); x1=0:471; y1=polyval(

7、a,x1); y2=polyval(b,x1); plot(x,y,*,x1,y1, r,x1,y2, g),血药浓度模型,通过实验测得一次性快速静脉注射300mg药物后的血药浓度数据见表：求血药浓度随时间的变化规律。,先画散点图,x=0.25 0.5 1 1.5 2 3 4 6 8; y=19.21 18.15 15.36 14.10 12.89 9.32 7.45 5.24 3.01; plot(x,y,*),根据散点图可知大致呈指数形态（y=a*exp(-bt)）,建立M文件（由于是单调递减，所以指数的系数为负的）function y=fun66(x,t) y=x(1)*exp(-x(2

8、)*t);%输入主程序 t=0.25 0.5 1 1.5 2 3 4 6 8; p=19.21 18.15 15.36 14.10 12.89 9.32 7.45 5.24 3.01; t0=10 0.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.2413 0.2420 p1 = Columns 1 through 8 19.0532 17.9348 15.8911 14.0802 12.4757 9.7945

9、7.6894 4.7394 Column 9 2.9211,建立M文件: function y=fun66(x,t) y=x(1)*exp(-x(2)*t);%输入主程序 t=0.25 0.5 1 1.5 2 3 4 6 8; p=19.21 18.15 15.36 14.10 12.89 9.32 7.45 5.24 3.01; t0=10 0.4292; x=nlinfit(t,p, fun66,t0) p1=fun66(x,t),拟合,t=0.25 0.5 1 1.5 2 3 4 6 8; p=19.21 18.15 15.36 14.10 12.89 9.32 7.45 5.24 3.

10、01; t1=0.25:0.01:8; p1=20.2413*exp(-0.2420*t1); plot(t,p,*) hold on plot(t1,p1,g:),t1=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:),转化为线性拟合来完成,t=0.25 0.5 1 1.5 2 3 4 6 8; y=19.21 18.15 15.36 14.10 12.89 9.32 7.45 5.24 3.01; Y=log(y); p=polyfit(t,Y,1) p = -0.2347 2.9943,即,syms t t=0.25:0.01:8; y=19.9714*exp(-0.2347*t); plot(t,y,g:) t0=0.25 0.5 1 1.5 2 3 4 6 8; p=19.21 18.15 15.36 14.10 12.89 9.32 7.45 5.24 3.01; hold on plot(t0,p,*),饮酒驾车的拟合图,t=0.25 0.5 0.75 1 1.5 2 2.5 3 3.5 4 4.5 5 6 7 8 9 10 11 12 13 14 15