【实验目的及要求】

熟练掌握dsolve，ode45，ode15s语句的调用格式；正确区分符号运算与数值运算的不同，避免出现混淆调用的情况；对实际问题，能够建立微分方程数学模型，能够用dsolve，ode45，ode15s等语句进行求解

【实验过程】(实验步骤、绘图、记录、数据、分析、结果及实验教师评语(可选))

一．实验题目

1.(必做题)解微分方程(组)

(1)

(提示可以考虑，以特解函数及其一阶、二阶导数曲线图形来表示)

解：①将高阶微分方程化为一阶微分方程组，设，则有

②建立函数文件

functiondy=myfun(x,y)

dy=[y(2);y(3);(y(3)-1)^2-y(2)-y(1)^2];

③主程序：

[x,y]=ode45('myfun',[0,20],[0;1;-1]);

plot(x,y(:,1),'*',x,y(:,2),'+',x,y(:,3),'o')

%legend('y','y的一阶导数','y的二阶导数');

④结果

注意此题得不到解析解，只能用数值解，解法可参看PPT中数值解例题3

(2)运用数值解手段描述下面常微分方程组在初值下的相空间的相轨线.

解：①建立函数文件

functiondx=lorenz(t,x)

dx=[-8/3*x(1)+x(2)*x(3);-10*x(2)+10*x(3);-x(1)*x(2)+28*x(2)-x(3)];

②主程序文件

[t,x]=ode45('lorenz',[0,100],[0;0;1e-10]);

axis([0 40 -20 20 -20 20]);

plot3(x(:,1),x(:,2),x(:,3));

gridon

③结果

(3)求的数值解，并画出图像

解：首先建立odefun1.m如下：

function dy=odefun1(x,y);

dy=[-0.01*y(1)-99.99*y(2);-100*y(2)];

然后建立主程序shiyan2_3.m

clc

clear

close all

[x,y]=ode15s('odefun1',[0 100],[2;1])

plot(x,y(:,1),'*',x,y(:,2),'r*')

结果为

(4)求下列方程的通解及特解

(Bessel方程，令n=½)

解：求通解的主程序为(syms n)

diff_y='x^2*D2y+x*Dy+(x^2-(1/2)^2)*y=0';

y=dsolve(diff_y,'x')

结果为：

ｙ＝C1/x^(1/2)*sin(x)+C2/x^(1/2)*cos(x)

y =

(2^(1/2)*C12*cos(x))/(pi^(1/2)*x^(1/2)) + (2^(1/2)*C13*sin(x))/(pi^(1/2)*x^(1/2))

求特解的主程序为

diff_y='x^2*D2y+x*Dy+(x^2-(1/2)^2)*y=0';

y=dsolve(diff_y,'y(pi/2)=2,Dy(pi/2)=-2/pi','x')

结果为：

y =2^(1/2)*pi^(1/2)/x^(1/2)*sin(x)

y =

(2*sin(x)*(pi/2)^(1/2))/x^(1/2) + (cos(x)*(2/(pi/2)^(1/2) - pi/(pi/2)^(3/2)))/(2*x^(1/2))

2.(必做题)凶杀案作案时间问题：受害者的尸体于晚上7:30被发现，法医于晚上8:20赶到凶案现场，测得尸体温度为32.6℃；一小时后，当尸体即将被抬走时，测得尸体温度为31.4℃，室温在几个小时内始终保持21.1℃。此案最大的嫌疑犯张某声称自己是无罪的，并有证人说：“下午张某一直在办公室上班，5:00时打完电话后就离开了办公室”。从张某到受害者家(凶案现场)步行需5分钟，现在的问题是，张某不在凶案现场的证言能否被采信，使他排除在嫌疑犯之外。(提示：Newton冷却定理告诉我们“物体在介质中冷却速度同该物体温度与介质温度之差成正比”)

解：首先应确定凶案的发生时间，若死亡时间在下午5点5分之前，则张某就不是嫌疑犯，否则不能将张某排除。

设T(t)表示t时刻尸体的温度，并记晚上8:20为t=0，则T(0)=32.6℃，T(1)=31.4℃。假设受害者死亡时体温是正常的，即T=37℃(查资料)是要确定受害者死亡的时间，也就是求T(t)=37℃的时刻，进而确定张某是否是嫌疑犯。

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

模型：由Newton冷却定理可得一阶线性微分方程模型

求解：(1)首先用dsolve求解该方程的解析解

程序：

symslamd

sy3d11='DT+lamd*(T-21.1)=0';

T=dsolve(sy3d11,'T(0)=32.6','t')

结果：

T =211/10+23/2*exp(-lamd*t)

或

T =23/(2*exp(lamd*t)) + 211/10

(2)求解参数lamd

可以利用初始条件“1小时后，当尸体即将被抬走时，测得尸体温度为31.4℃”

由上式可以得到：31.4=21.1+11.5*exp(-1*lamd)

lamd的值为0.11020314013361429463890984998294

(程序：lamd=solve('31.4-21.1-11.5*exp(-1*lamd)=0','lamd')

(3)求解t0

当T=37℃时，有t=-2.95小时

(程序：t0=solve('37-21.1-11.5*exp(-0.11*t)','t'))

＝-2小时57分，

8小时20分－2小时57分＝5小时23分。即死亡时间大约在5:23，

因此张某不能被排除在嫌疑犯之外。