1、数学数学上机教学（三）上机教学（三）微分方程求解微分方程求解上机目的上机目的上机内容上机内容、学会用、学会用 求微分方程的数值解求微分方程的数值解.上机软件上机软件、学会用、学会用 求简单微分方程的解析解求简单微分方程的解析解.、求简单微分方程的解析解、求简单微分方程的解析解.、上机作业、上机作业.、求微分方程的数值解、求微分方程的数值解.、数学建模实例数学建模实例.、微分方程的解析解、微分方程的解析解 求微分方程（组）的解析解命令求微分方程（组）的解析解命令:dsolve(方程方程1,方程方程2,方程方程n,初始条件初始条件,自变量自变量)结果：结果：()例如求下例微分方程的特解。例如求下例
2、微分方程的特解。,(0)xdyeyedx在命令窗口中输入：在命令窗口中输入：(),()(),)输出结果为：输出结果为：()()如想画出函数在自变量在区间如想画出函数在自变量在区间的函数图像，可在的函数图像，可在命令窗口中输入：命令窗口中输入：(,)解解 输入命令输入命令:(*,()(),)结结 果果 为为:*(*)*(*)作图命令：作图命令：(,)解解 输入命令输入命令：x,y,z=dsolve(Dx=2*x-3*y+3*z,Dy=4*x-5*y+3*z,Dz=4*x-4*y+2*z,t);x=simple(x)%将将x化简化简 y=simple(y)z=simple(z)运行结果为：运行结果
3、为：x=(c1-c2+c3+c2e-3t-c3e-3t)e2t,y=-c1e-4t+c2e-4t+c2e-3t-c3e-3t+c1-c2+c3)e2t ,z=(-c1e-4t+c2e-4t+c1-c2+c3)e2t.、微分方程的数值解、微分方程的数值解（一）常微分方程数值解的定义（一）常微分方程数值解的定义 在生产和科研中所处理的微分方程往往很复杂且大多得在生产和科研中所处理的微分方程往往很复杂且大多得不出一般解不出一般解.而在实际上对初值问题，一般是要求得到解在若而在实际上对初值问题，一般是要求得到解在若干个点上满足规定精确度的近似值，或者得到一个满足精确干个点上满足规定精确度的近似值，或者
4、得到一个满足精确度要求的便于计算的表达式度要求的便于计算的表达式.因此，研究常微分方程的数值解法是十分必要的因此，研究常微分方程的数值解法是十分必要的.0000121212yf(x,y)x y(x)y(x),(),()y,nnnyxxxxxy xy xyy对常微分方程：，其数值解是指由初始点开始的若干离散的 值处，即对，求出准确值的相应近似值。（二）建立数值解法的一些途径（二）建立数值解法的一些途径001i)y(xy)f(x,y ,1,2,1,0 ,xynihxi解微分方程：可用以下离散化方法求设、用差商代替导数、用差商代替导数 若步长若步长h较小，则有较小，则有hxyhxyxy)()()(故
5、有公式：故有公式：100(,),i0,1,2,n-1.(),iiiiyyhf xyyy x此即欧拉法此即欧拉法.、使用数值积分、使用数值积分对方程对方程y=f(x,y),两边由两边由xi到到xi+1积分，并利用梯形公式，有：积分，并利用梯形公式，有：)(,()(,(2)(,()()(11111iiiiiixxiixyxfxyxfxxdttytfxyxyii实际应用时，与欧拉公式结合使用：实际应用时，与欧拉公式结合使用：,2,1,0 ),(),(2),()(11)1(1)0(1kyxfyxfhyyyxhfyykiiiiikiiiii的计算。然后继续下一步，取时，当满足，对于已给的精确度）（y y
6、 2i111i)(1)1(1kikikiyyy此即改进的欧拉法此即改进的欧拉法.故有公式：故有公式：)(),(),(200111xyyyxfyxfhyyiiiiii、使用泰勒公式、使用泰勒公式 以此方法为基础，有龙格库塔法、线性多步法等方法以此方法为基础，有龙格库塔法、线性多步法等方法.、数值公式的精度、数值公式的精度 当一个数值公式的截断误差可表示为（）时（为正整当一个数值公式的截断误差可表示为（）时（为正整数，为步长），称它是一个阶公式数，为步长），称它是一个阶公式.越大，则数值公式的精度越高越大，则数值公式的精度越高.欧拉法是一阶公式，改进的欧拉法是二阶公式欧拉法是一阶公式，改进的欧拉法
7、是二阶公式.龙格库塔法有二阶公式和四阶公式龙格库塔法有二阶公式和四阶公式.线性多步法有四阶阿达姆斯外插公式和内插公式线性多步法有四阶阿达姆斯外插公式和内插公式.（三）用软件求常微分方程的数值解（三）用软件求常微分方程的数值解，（）ode45 ode23 ode113ode15sode23s由待解由待解方程写方程写成的成的m-文件名文件名ts=t0，tf，t0、tf为自为自变量的初变量的初值和终值值和终值函数的函数的初值初值ode23：组合的：组合的2/3阶龙格阶龙格-库塔库塔-芬尔格算法芬尔格算法ode45：运用组合的：运用组合的4/5阶龙格阶龙格-库塔库塔-芬尔格算法芬尔格算法自变自变量值量
8、值函数函数值值用于设定误差限用于设定误差限(缺省时设定相对误差缺省时设定相对误差10-3,绝对误差绝对误差10-6),命令为：命令为：options=odeset（reltol,rt,abstol,at）,rt，at：分别为设定的相对误差和绝对误差：分别为设定的相对误差和绝对误差.、在解个未知函数的方程组时，和均为维向量，文件中的、在解个未知函数的方程组时，和均为维向量，文件中的待解方程组应以待解方程组应以 的分量形式写成的分量形式写成.、使用、使用 软件求数值解时，高阶微分方程必须等价地变换软件求数值解时，高阶微分方程必须等价地变换成一阶微分方程组成一阶微分方程组.注意注意:解解:令令 y1
9、=x，y2=y1、建立文件如下：、建立文件如下：()();()();()*()*()();、取，在主命令窗口中输入命令：、取，在主命令窗口中输入命令：(,);(),)、结果如图：、结果如图：050010001500200025003000-2.5-2-1.5-1-0.500.511.52解解 、建立文件如下：、建立文件如下：()();()()*();()()*();()*()*();、取，输入命令：、取，输入命令：(,);(),(),*(),)、结果如图：、结果如图：024681012-1-0.8-0.6-0.4-0.200.20.40.60.81 图中，的图形为实线，的图形为图中，的图形为实
10、线，的图形为“*”线，的图形为线，的图形为“”“”线线.（一）导弹追踪问题（一）导弹追踪问题 设位于坐标原点的甲舰向位于设位于坐标原点的甲舰向位于x轴上点轴上点A(1,0)处的乙舰处的乙舰发射导弹，导弹头始终对准乙舰发射导弹，导弹头始终对准乙舰.如果乙舰以最大的速度如果乙舰以最大的速度v0(是常数是常数)沿平行于沿平行于y轴的直线行驶，导弹的速度是轴的直线行驶，导弹的速度是5v0，求，求导弹运行的曲线方程导弹运行的曲线方程.又乙舰行驶多远时，导弹将它击中？又乙舰行驶多远时，导弹将它击中？解法一：（解析法）解法一：（解析法）假设导弹在t 时刻的位置为P(x(t),y(t)，乙舰位于),1(0tv
11、Q.由于导弹头始终对准乙舰，故此时直线 PQ就是导弹的轨迹曲线弧 OP 在点 P 处的切线，即有 xytvy10即 yyxtv)1(0 (1)又根据题意，弧 OP 的长度为AQ的 5 倍，即 tvdxyx00251 (2)3、数学建模实例、数学建模实例由由(1),(2)消去消去t整理得模型整理得模型:21(1)1,(3)5x yy解即为导弹的运行轨迹:245)1(125)1(855654xxy当1x时245y，即当乙舰航行到点)245 ,1(处时被导弹击中.被击中时间为:00245vvyt.若 v0=1,则在 t=0.21 处被击中.解法二：解法二：(数值解数值解).建立文件建立文件 ()()
12、;()();()*()();.取取;，建立主程序如下：，建立主程序如下：;(,);(),.);(,*)结论结论:导弹大致在（，）处击中乙舰导弹大致在（，）处击中乙舰.21(1)15x yy1222211/(1)5yyyyx令令，将方程（）化为一阶微分方程组，将方程（）化为一阶微分方程组.解法三：解法三：(建立参数方程求数值解建立参数方程求数值解)设时刻设时刻t乙舰的坐标为乙舰的坐标为(X(t)，Y(t)，导弹的坐标为，导弹的坐标为(x(t)，y(t).1设导弹速度恒为w，则 222)()(wdtdydtdx （1）2.由于弹头始终对准乙舰，故导弹的速度平行于乙舰与导弹头位置的差向量，即：yYx
13、Xdtdydtdx，0 （2）消去得：)()()()()()(2222yYyYxXwdtdyxXyYxXwdtdx （3）3因乙舰以速度因乙舰以速度v0沿直线沿直线x=1运动，设运动，设v0=1，则，则w=5，X=1，Y=t4.解导弹运动轨迹的参数方程解导弹运动轨迹的参数方程(1)建立建立m-文件文件eq2.m如下：如下：function dy=eq2(t,y)dy=zeros(2,1);dy(1)=5*(1-y(1)/sqrt(1-y(1)2+(t-y(2)2);dy(2)=5*(t-y(2)/sqrt(1-y(1)2+(t-y(2)2);（2）取）取t0=0，tf=2，建立主程序，建立主程
14、序chase2.m如下：如下：t,y=ode45(eq2,0 2,0 0);Y=0:0.01:2;plot(1,Y,-)hold on;plot(y(:,1),y(:,2),*)5.结果见图结果见图1导弹大致在（导弹大致在（1，0.2）处击中乙舰，与前面的结论一致）处击中乙舰，与前面的结论一致.图图1图图2 在在chase2.m中，按二分法逐步修改中，按二分法逐步修改tf，即分别取，即分别取tf=1，0.5,0.25,直到直到tf=0.21时，得图时，得图2.结论：时刻结论：时刻t=0.21时，导弹在（时，导弹在（1，0.21）处击中乙舰）处击中乙舰.、上、上 机机 作作 业业（三）（三）.鱼
15、雷追击问题鱼雷追击问题 一敌舰在某海域内沿着正北方向航行时，我方战舰恰好位于一敌舰在某海域内沿着正北方向航行时，我方战舰恰好位于敌舰的正西方向敌舰的正西方向 公里处我舰向敌舰发射制导鱼雷，敌舰速度为公里处我舰向敌舰发射制导鱼雷，敌舰速度为 公里分，鱼雷速度为敌舰速度的倍。试问敌舰航行多远时将被击公里分，鱼雷速度为敌舰速度的倍。试问敌舰航行多远时将被击中？中？.求微分方程求微分方程 ，在初值条件，在初值条件 下的下的特解，并画出解函数的图形特解，并画出解函数的图形.0 xxyye 12()ye 一个慢跑者在平面上沿椭圆以恒定的速率一个慢跑者在平面上沿椭圆以恒定的速率v=1跑跑步步,设椭圆方程为设椭圆方程为:x=10+20cost,y=20+5sint.突然有一突然有一只狗攻击他只狗攻击他.这只狗从原点出发这只狗从原点出发,以恒定速率以恒定速率w跑向慢跑者跑向慢跑者,狗的运动方向始终指向慢跑者狗的运动方向始终指向慢跑者.分别求出分别求出w=20,w=5时狗的时狗的运动轨迹运动轨迹.(选作选作)2009年6月（慢跑者与狗）（慢跑者与狗）
