实验4常微分方程数值解

《实验4常微分方程数值解》由会员上传分享，免费在线阅读，更多相关内容在教育资源-天天文库。

1、实验4常微分方程数值解凯里学院理学院主讲:潘东云ExperimentsinMathematics1为什么要学习微分方程数值解微分方程是研究函数变化规律的重要工具，有着广泛的应用。如:物体的运动,电路的电压,人口增长的预测许多微分方程没有解析解，数值解法是求解的重要手段，如2实验4的基本内容3.实际问题用微分方程建模，并求解2.龙格-库塔方法的MATLAB实现*4.数值算法的收敛性、稳定性与刚性方程两个最常用的数值算法:欧拉（Euler）方法龙格-库塔（Runge-Kutta）方法3实例1海上缉私海防某部缉私艇上的雷达发现正东方向c海里处有一艘走私船正以速度

2、a向正北方向行驶，缉私艇立即以最大速度b(>a)前往拦截。如果用雷达进行跟踪时，可保持缉私艇的速度方向始终指向走私船。建立任意时刻缉私艇位置及航线的数学模型,并求解;求出缉私艇追上走私船的时间。a北bc艇船4实例1海上缉私建立坐标系如图:t=0艇在(0,0),船在(c,0);船速a,艇速b时刻t艇位于P(x,y),船到达Q(c,at)模型:0yxcR(c,y)Q(c,at)P(x,y)b由方程无法得到x(t),y(t)的解析解需要用数值解法求解5“常微分方程初值问题数值解”的提法不求解析解y=y(x)(无解析解或求解困难)而在一系列离散点通常取等步长h

3、y1y2yn6yP0x0x1x2x3xy=y(x)y0P1P2P3欧拉方法基本思路x取不同点各种欧拉公式向前欧拉公式显式公式7P3欧拉方法向后欧拉公式隐式公式yP0x0x1x2x3xy0y=y(x)P1P28向前欧拉公式向后欧拉公式二者平均得到梯形公式仍为隐式公式，需迭代求解将梯形公式的迭代过程简化为两步预测校正改进欧拉公式9龙格-库塔方法向前，向后欧拉公式：用[xn,xn+1]内１个点的导数代替f(x,y(x))梯形公式，改进欧拉公式：用[xn,xn+1]内２个点导数的平均值代替f(x,y(x))龙格-库塔方法的基本思想在[xn,xn+1]内多取几个点，

4、将它们的导数加权平均代替f(x,y(x))，设法构造出精度更高的计算公式。10常用的(经典)龙格—库塔公式不足:收敛速度较慢L级?阶龙格-库塔方法的一般形式使精度尽量高4级4阶11微分方程组和高阶方程初值问题的数值解欧拉方法和龙格-库塔方法可直接推广到微分方程组向前欧拉公式高阶方程需要先降阶为一阶微分方程组12龙格—库塔方法的MATLAB实现[t,x]=ode23(@f,ts,x0,opt)3级2阶龙格-库塔公式[t,x]=ode45(@f,ts,x0,opt)5级4阶龙格-库塔公式f是待解方程写成的函数m文件：functiondx=f(t,x)dx=[f

5、1;f2;;fn];ts=[t0,t1,…,tf]输出指定时刻t0,t1,…,tf的函数值ts=t0:k:tf输出[t0,tf]内等分点处的函数值x0为函数初值(n维)输出t=ts,x为相应函数值(n维)opt为选项，缺省时精度为：相对误差10-3，绝对误差10-6,计算步长按精度要求自动调整13微分方程的符号解求微分方程（组）的解析解(符号解）命令:dsolve(‘方程1’,‘方程2’,…‘方程n’,‘初始条件’,‘自变量’)注意:在表达微分方程时，用字母D表示求微分，D2、D3等表示求高阶微分.任何D后所跟的字母为因变量，自变量可以指定或由系统规则选

6、定为确省.例如，微分方程应表达为：D2y=0.结果：u=tan(t-c)14解输入命令:y=dsolve('D2y+4*Dy+29*y=0','y(0)=0,Dy(0)=15','x')结果为:y=3e-2xsin（5x）15解输入命令：[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)结果为：x=(c1-c2+c3+c2e-3t-c3e-3t)e2ty=-c1e-4t+c2e-4t+c2e

7、-3t-c3e-3t+c1-c2+c3)e2tz=(-c1e-4t+c2e-4t+c1-c2+c3)e2t16解：用MATLAB编程如下：functiondx=fun113(t,x)dx=t+x;ts=0:0.1:1;x0=1;[t,x]=ode45('fun113',ts,x0);y=2*exp(t)-t-1;[t,x,y,x-y]plot(t,x,'r-',t,y,'b-.'),grid,gtext('x(t)'),gtext('y(t)'),title('微分方程数值解例');xlabel('timet');ylabel('solutionx');l

8、egend('x','y');微分方程的数值解17微分方程的数值解