《数学实验 第4版》课件 4.4 常微分方程的数值解

第四章数值分析

实验4.1插值

实验4.2离散数据的曲线拟合数学实验实验4.3MATLAB数值积分与微分实验4.4常微分方程的数值解2实验4.4常微分方程的数值解一、几种求常微分方程数值解的方法二、应用举例实验4.4常微分方程的数值解3实验4.4常微分方程的数值解实验目的通过本实验了解常微分方程的数值解的概念，掌握利用MATLAB求常微分方程的数值解的方法.一、几种求常微分方程数值解的方法常微分方程是研究函数变化规律的有力工具，但是绝大多数变系数方程、非线性方程都是所谓“解不出来”的，于是常微分方程的数值解法就成为解常微分方程的主要手段．本实验只考虑初值问题.4实验4.4常微分方程的数值解常微分方程初值问题：保证方程的解存在且惟一在一系列离散点上，求的近似值，通常取等步长h，即因此数值解法得到的近似解是一个离散的函数表.5实验4.4常微分方程的数值解一、几种求常微分方程数值解的方法1.欧拉（Euler）方法欧拉方法是一种最古老、最简单而直观的解微分方程的数值方法，其基本想法是在小区间上用差商代替方程左端的导数,而方程右端已知函数中的x在小区间上的哪一点取值，则有以下不同的方法：6实验4.4常微分方程的数值解（1）向前欧拉公式中的x取小区间的左端点xn，以分别代替，就得到：向前欧拉公式（2）向后欧拉公式中的x取小区间的右端点xn+1，就得到：向后欧拉公式上式右端的yn+1未知，故称为隐式公式，无法用它直接计算yn+1．7实验4.4常微分方程的数值解（3）梯形公式将向前欧拉公式和向后欧拉公式加以平均，得到梯形公式8实验4.4常微分方程的数值解（4）改进的欧拉公式先由向前欧拉公式算出yn+1的预测值，再把它代入梯形公式右端，作为校正，即改进的欧拉公式9实验4.4常微分方程的数值解它还可写作上面4个公式中，我们常用的是便于计算的向前欧拉公式和改进的欧拉公式.10实验4.4常微分方程的数值解2.龙格—库塔（Runge-Kutta）公式在区间内多取几个点，就可以构造出精度更高的计算公式，这就是龙格—库塔公式，它是一类方法的总称.它是欧拉方法的一种推广，也是应用最广的求解常微分方程数值问题的方法.一般的龙格—库塔方法的形式为p阶龙格—库塔方法其中ai,bij,ci为待定参数11实验4.4常微分方程的数值解（1）二阶龙格—库塔公式其中为待定系数.满足：上式有4个未知数而只有3个方程，所以解不唯一，存在无穷多个解，可见二阶龙格—库塔方法是一族公式.

改进的欧拉公式12实验4.4常微分方程的数值解（2）四阶经典（标准）龙格—库塔公式实际应用中，用的最多的是四阶龙格—库塔公式四阶龙格—库塔公式也不止一个13实验4.4常微分方程的数值解3.龙格—库塔方法的MATLAB实现在MATLAB中，求微分方程数值解的函数调用格式为[t,x]=solver（’f’,ts,x0,options）自变量值函数值ode45、ode23、ode113、ode15s、ode23s中的一个由待解方程写成的m-文件名=[t0,tf],t0、tf为自变量的初值和终值函数的初值用于设定误差限options=odeset（'reltol',rt,'abstol',at）默认相对误差为10-3，绝对误差为10-6.相对误差绝对误差14实验4.4常微分方程的数值解在诸多的solver函数中，ode23、ode45都是基于龙格—库塔公式的函数:ode23是采用的2、3阶龙格-库塔组合算法，ode45是采用的4、5阶龙格-库塔组合算法.ode45比ode23精度更高，尤其成为大部分场合的首选算法.注(1)在解n个未知函数的方程组时，x0和x均为n维向量，m-文件中的待解方程组应以x的分量形式写成.(2)使用Matlab软件求数值解时，高阶微分方程必须等价地变换成一阶微分方程组.15例21

求微分方程的数值解，解先编写函数文件：functiondx=fun23(t,x)dx=t+x;再求方程的解实验4.4常微分方程的数值解已知精确解为>>ts=0:0.1:1;x0=1;[t,x]=ode45('fun23',ts,x0);y=2*exp(t)-t-1;[t,x,y]plot(t,x,'r-',t,y,'b-.'),gridon;title('SolutionofExample3');xlabel('timet');ylabel('solutionx');legend('x','y');↙ans=01.00001.000016实验4.4常微分方程的数值解0.10001.11031.11030.20001.24281.24280.30001.39971.39970.40001.58361.58360.50001.79741.79740.60002.04422.04420.70002.32752.32750.80002.65112.65110.90003.01923.01921.00003.43663.4366

t的值x的值y的值例22

求微分方程组的数值解.实验4.4常微分方程的数值解17解先编写函数文件：functiondx=fun24(t,x)dx=[3*t*x(1)+6*x(2);t^2+x(1)*x(2)];再求方程组的解>>ts=0:0.1:1.5;x0=[1,0];[t,x]=ode45('fun24',ts,x0);a=[t,x];plot(t,x(:,1),'r-',t,x(:,2),'b-.'),gridon↙18a=01.000000.10001.01520.00030.20001.06270.00280.30001.14910.00980.40001.28650.02410.50001.49520.04940.60001.80790.09110.70002.27750.15850.80002.99200.27010.90004.10490.4707实验4.4常微分方程的数值解1.00005.90410.88741.10009.01972.00691.200015.32486.75821.300037.960170.5533t的值x的值y的值实验4.4常微分方程的数值解例23求解描述振荡器的经典的VerderPol微分方程:解原方程等价为，，.，1920首先编写函数文件verderpol.m：functionxprime=verderpol(t,x)globalmu;xprime=[x(2);mu*(1-x(1)^2)*x(2)-x(1)]；再在命令窗口中执行：globalmu;mu=7;y0=[1;0];[t,x]=ode45('verderpol',[0,40],y0);x1=x(:,1);x2=x(:,2);plot(t,x1,t,x2)↙实验4.4常微分方程的数值解21通过这些繁杂的公式，我们看到数学家们为了模型的完善和改进付出无数的心血和努力，这种对科学精益求精的探索和创新精神值得我们学习.同时，为了更好地应用这些成果，我们应该具体问题具体分析，选取适当的方法解决问题.实验4.4常微分方程的数值解22实验4.4常微分方程的数值解二、应用举例导弹追踪问题设位于坐标原点的甲舰向位于x轴上点A(1,0)处的乙舰发射导弹，导弹头始终对准乙舰．如果乙舰以最大的速度v0(是常数)沿平行于y轴的直线行驶，导弹的速度是5v0，求导弹运行的曲线方程．又乙舰行驶多远时，导弹将它击中？解法一（解析法）：假设导弹在t时刻的位置为乙舰位于.由于导弹头始终对准乙舰，故此时直线PQ就是导弹的轨迹曲线弧在点P处的切线，如图所示：23实验4.4常微分方程的数值解图4.14即有即(1)又根据题意,弧的长度为的5倍，即(2)由(1)，(2)两式消去t整理得模型:初值条件为:初值问题的解即为导弹的运行轨迹:24当时,即当乙舰航行到点处时被导弹击中.实验4.4常微分方程的数值解若,则在处被击中.被击中时间为:编写以下程序可得轨迹图>>x=0:0.01:1;y=-5*(1-x).^(4/5)/8+5*(1-x).^(6/5)/12+5/24;plot(x,y,'*')↙25实验4.4常微分方程的数值解解法二(数值解)(1)建立m-文件eq1.mfunctiondy=eq1(x,y)dy=zeros(2,1);dy(1)=y(2);dy(2)=1/5*sqrt(1+y(1)^2)/(1-x);(2)取x0=0，xf=0.9999，求解如下：>>x0=0;xf=0.9999;[x,y]=ode15s('eq1',[x0xf],[00]);plot(x,y(:,1),'b.')holdony=0:0.01:2;plot(1,y,'b*')↙26结论:导弹大致在（1,0.2）处击中乙舰.实验4.4常微分方程的数值解27实验4.4常微分方程的数值解解法三(建立参数方程求数值解)设时刻乙舰的坐标为，导弹的坐标为(1)设导弹速度恒为(2)由于弹头始终对准乙舰，故导弹的速度平行于乙舰与导弹头位置的差向量,即：28实验4.4常微分方程的数值解消去λ得：(3)因乙舰以速度沿直线运动因此导弹运动轨迹的参数方程为：，则设29实验4.4常微分方程的数值解(4)解导弹运动轨迹的参数方程建立m-文件eq2.m如下：functiondy=eq2(t,y)dy=zeros(2,1);dy(1)=5*(1-y(1))/sqrt((1-y(1)