大学数学实验.ppt

《大学数学实验.ppt》由会员上传分享，免费在线阅读，更多相关内容在工程资料-天天文库。

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

2、东方向c海里处有一艘走私船正以速度a向正北方向行驶，缉私艇立即以最大速度b(>a)前往拦截。如果用雷达进行跟踪时，可保持缉私艇的速度方向始终指向走私船。建立任意时刻缉私艇位置及航线的数学模型,并求解;求出缉私艇追上走私船的时间。a北bc艇船实例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)的解析解需要用数值解法求解“常微分方程初值问题数值解”的提法不求解析解y=y(x)(无

3、解析解或求解困难)而在一系列离散点通常取等步长hy1y2ynyP0x0x1x2x3xy=y(x)y0P1P2P3欧拉方法基本思路x取不同点各种欧拉公式向前欧拉公式显式公式P3欧拉方法向后欧拉公式隐式公式yP0x0x1x2x3xy0y=y(x)P1P2向前欧拉公式向后欧拉公式二者平均得到梯形公式仍为隐式公式，需迭代求解将梯形公式的迭代过程简化为两步预测校正改进欧拉公式龙格-库塔方法向前，向后欧拉公式：用[xn,xn+1]内１个点的导数代替f(x,y(x))梯形公式，改进欧拉公式：用[xn,xn+1]内２个点导数的平均值代替f(x,y(x

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

5、数m文件：functiondx=f(t,x)dx=[f1;f2;;fn];ts=[t0,t1,…,tf]输出指定时刻t0,t1,…,tf的函数值ts=t0:k:tf输出[t0,tf]内等分点处的函数值x0为函数初值(n维)输出t=ts,x为相应函数值(n维)缺省精度(相对误差10-3，绝对误差10-6),计算步长按精度要求自动调整.实例1海上缉私(续)模型的数值解0yxc(x(t),y(t))ab设:船速a=20(海里/小时)艇速b=40(海里/小时)距离c=15(海里)求:缉私艇的位置x(t),y(t)缉私艇的航线y(x)实例1海上

6、缉私(续)模型的数值解a=20,b=40,c=15走私船的位置x1(t)=c=15y1(t)=at=20tt=0.5时缉私艇追上走私船缉私艇的航线y(x)tx(t)y(t)0000.051.99840.06980.103.98540.29240.155.94450.69060.207.85151.28990.259.67052.11780.3011.34963.20050.3512.81704.55520.4013.98066.17730.4514.74518.02730.5015.00469.9979y1(t)01.02.03.04.

7、05.06.07.08.09.010.0刚性现象与刚性方程振动系统或电路系统的数学模型现象k=2000.5,r=1000,a=1,b=-1999.5,f(t)=1瞬态解与稳态解e-2000t~快瞬态解e-t/2~慢瞬态解计算到t=0.005时已衰减到4.510-5计算到t=20时才衰减到4.510-5精度达到10-4需算到t=20(由慢瞬态解=1/2决定)选取步长h由快瞬态解=2000决定h<2.785/2000=0.0014龙格-库塔公式t=20需14286步快、慢瞬态解的特征根相差悬殊刚性现象(Stiff)求稳态解刚性现象与

8、刚性方程刚性方程的MATLAB求解ode23,ode45解刚性方程的困难步长自动变小计算时间很长求解刚性方程的命令:ode23s,ode15s等(用法相同)例特征根1=-1,2=-106刚性比s=106