专题四 微分方程的matlab求解
2012-12-15 20:15阅读：
【实验目的及要求】
熟练掌握dsolve，ode45，ode15s语句的调用格式；正确区分符号运算与数值运算的不同，避免出现混淆调用的情况；对实际问题，能够建立微分方程数学模型，能够用dsolve，ode45，ode15s等语句进行求解
【实验过程】（实验步骤、绘图、记录、数据、分析、结果及实验教师评语（可选））
一．实验题目
1.（必做题）解微分方程（组）
（1）
![专题四 <wbr>微分方程的matlab求解]( "专题四 <wbr>微分方程的matlab求解")
（提示可以考虑 ![专题四 <wbr>]()
微分方程的matlab求解' TITLE='专题四
微分方程的matlab求解'
/>
，以特解函数及其一阶、二阶导数曲线图形来表示）
解：
①将高阶微分方程化为一阶微分方程组，设
![专题四 <wbr>微分方程的matlab求解]( "专题四 <wbr>微分方程的matlab求解")
，则有
![专题四 <wbr>微分方程的matlab求解]( "专题四 <wbr>微分方程的matlab求解")
②建立函数文件
function
dy=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的二阶导数');
④结果
![专题四 <wbr>微分方程的matlab求解]( "专题四 <wbr>微分方程的matlab求解")
注意此题得不到解析解，只能用数值解，解法可参看PPT中数值解例题3
（2）运用数值解手段描述下面常微分方程组在初值
![专题四 <wbr>微分方程的matlab求解]( "专题四 <wbr>微分方程的matlab求解")
下的相空间的相轨线.
![专题四 <wbr>微分方程的matlab求解]( "专题四 <wbr>微分方程的matlab求解")
解：①建立函数文件
function dx=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));
grid on
③结果
![专题四 <wbr>微分方程的matlab求解]( "专题四 <wbr>微分方程的matlab求解")
（3）求 ![专题四 <wbr>微分方程的matlab求解]( "专题四 <wbr>微分方程的matlab求解")
的数值解，并画出图像
解：首先建立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*')
结果为 ![专题四 <wbr>微分方程的matlab求解]( "专题四 <wbr>微分方程的matlab求解")
（4）求下列方程的通解及特解
![专题四 <wbr>微分方程的matlab求解]( "专题四 <wbr>微分方程的matlab求解")
（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求解该方程的解析解
程序：
syms lamd
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，
因此张某不能被排除在嫌疑犯之外。
