计算机仿真运算法

数学实验Experiments in Mathematics Laboratory Mathematics阮小娥博士李萍

2014/3/31

计算机仿真运算法

数学实验报告要求:1.实验问题 2.问题分析:包括解决问题的理论依据,建立的数学模型以及求解问题的思路和方法. 3.程序设计. 4.结果分析和结论. 5.总结和体会.

计算机仿真运算法

实验3缉私艇追击走私船3.1实验目的

2014/3/31

(1)学会用MATLAB软件求解微分方程的初值问题; (2)了解微分方程数值解的思想,掌握微分方程数值解的方法; (3)学会根据实际问题建立简单微分方程数学模型,提高解决问题能力; (4)了解简单的计算机仿真和数据模拟的基本方法. 3.2理论知识 (1)微分方程及解析解; (2)微分方程及数值解;在实际生产和科研中所建立的微分方程往往很复杂,尽管是一阶微分方程,有时也很难求解,有些微分方程甚至只能求出其隐式解，要想计算相应已知变量的函数值也是很困难的,在这种情况下，我们一般转向去求微分方程及数值解，即近似解。在实际上，对初值问题，一般是要求得到解在若干个点上满足规定精确度的近似值，研究常微分方程的数值解法，或者得到一个满足精确度要求的便于计算的表达式。

计算机仿真运算法

因此，研究常微分方程的数值解法是十分必要的。对常微分方程 y f ( x, y ) y ( x0 ) y0其数值解就是求其解析解

y y ( x)

在一系列离散点

x1 x2 xn xn 1 处的近似值

y1, y 2, , y n,

计算机仿真运算法

3.3实验问题 1用MATLAB软件求解析解MATLAB软件提供的解常微分方程解析解的指令是dsolve,完整的调用格式是:

dsolve(‘方程1’,‘方程2’,…‘方程n’,‘初始条件’,‘自变量’)微分方程的书写格式规定：当y是因变量时，用“Dny”表示y的n阶导数。

例1求微分方程

y x xy

的通解。

解

输入命令: y=dsolve('Dy=x+x*y','x') y=-1+exp(1/2*x^2)*C1

计算机仿真运算法

例2

求微分方程的特解.

d 2 y dy 2 4 29 y 0 dx dx y (0) 0, y ' (0) 15解输入命令: y=dsolve('D2y+4*Dy+29*y=0','y(0)=0,Dy(0)=15','x')

结果为: y=3e-2xsin(5x)

计算机仿真运算法

2用MATLAB软件求微分方程的数值解一般格式为

[t,x]=solver(’f’,ts,x0,options)自变量值函数值 ode45 ode23 ode113 ode15s ode23s由待解方程写成的m文件名 ts=[t0, tf],t0、 tf为自变量的初值和终值函数的初值

ode23:组合的2/3阶龙格-库塔-芬尔格算法 ode45:运用组合的4/5阶龙格-库塔-芬尔格算法用于设定误差限(缺省时设定相对误差10-3,绝对误差10-6,命令为：options=odeset(’reltol’,rt,’abstol’,at), rt,at:分别为设定的相对误差和绝对误差.注意:1、在解n个未知函数的方程组时，x0和x均为n维向量，m-文件中的待解方程组应以x的分量形式写成. 2、使用Matlab软

件求数值解时，高阶微分方程必须等价地变换成一阶微分方程组.

计算机仿真运算法

例1 y x, y (0) 0 function u=文件名(自变量,因变量) u=微分方程的右端然后在指令窗口执行 ode23(‘文件名’,[0,2],0)建立M文件,格式如左

可得到图形与相关数据可得到两组数据

再执行[x,y]=ode23(‘文件名’,[0,2],0)执行plot(x,y)可得到图形

polyfit(自变量，因变量，次数)

计算机仿真运算法

例2

y 2 y x 2 y (0) 1(1)求解析解; (2)求数值解,画出数值解与解析解的曲线比较.实验过程:

(1) y=dsolve('Dy=2*y+x+2','y(0)=1','x')结果为y=-1/2*x-5/4+9/4*exp(2*x) x=0:0.05:6; y=-1/2*x-5/4+9/4*exp(2*x); plot(x,y)

计算机仿真运算法

(2)编写函数文件,并存为fc10 function f=fc10(x,y) f=2*y+x+2然后在指令窗口执行[xb,yb]=ode23('fc10',[0,1],1) plot(x,y,xb,yb,’r*’)

计算机仿真运算法

3缉私艇追击走私船海上边防缉私艇发现距c公里处有一走私船正以匀速a沿直线行驶,缉私艇立即以最大速度b追赶,在雷达的引导下，缉私艇的方向始终指向走私船。问缉私艇何时追赶上走私船?并求出缉私艇追赶的路线。

y

o

c

x

计算机仿真运算法

建立模型走私船初始位置在点(0,0),行驶方向为y轴正方向，缉私艇的初始位置在点(c,0),缉私艇行驶的路程为s。在时刻t:走私船的位置到达点 R(0, at )缉私艇到达点 D ( x, y )

(0, at )

( x, y )

o2 d2y dy x 2 r 1 dx dx y (c) 0, y (c) 02

dy y at tg dx x 0

ds b dt

d y dt x a 2 dx dx

2

dt dt ds 1 dy 1 dx ds dx b dx

r a/b

计算机仿真运算法

模型求解(1)求解析解

d 2 y dp dy p, 2 令： dx, dx dx dx dp r 2 x 1 p p (c ) 0 r x 2 p 1 p c r c p 1 p 2 x

2 d2y dy x 2 r 1 dx dx y (c ) 0, y (c ) 0

r a/b dy 1 x r c r dx 2 c x y (c ) 0

计算机仿真运算法

a 1)r 1, b

1 r 1 r c 1 x 1 x cr y 2 2 1 r c 1 r c 1 r y cr bc cr 2当 x= 0时， y 2, t 2 2 1 r a a (1 r ) (b a )

dy 1 x r c r dx 2 c x y (c ) 0

c=3千米，a=0.4千米/分，分别取b=0.6,0.8,1.2千米/分时，缉私艇追赶路线的图形。

4 3.5 3 2.5 2

追赶时间分别为： t=9,5,2.8125(分钟)

1.5 1 0.5 0 0 0.5 1 1.5 2 2.5 3 3.5

计算机仿真运算法

a 2) r 1, b

1 r r 1 1 c c 1 x cr y 2 2 r 1 x r 1 1 r c 当 x 0

时， y ,缉私艇不可能追赶上走私船。

dy 1 x r c r dx 2 c x y (c ) 0

3) r

1

1 x2 c2 x , y 2 2c c ln c y ,缉私艇不可能追赶上走私船。

当 x 0时，

计算机仿真运算法

(2)用MATLAB软件求数值解c=3,a=0.4,b=0.8, r a/ b 0.5程序zjwt.m function u=zjwt(t,y) u=0.5*((t/3)^0.5-(3/t)^0.5)执行下面的命令：

dy 1 x r c r dx 2 c x y (c ) 0

ode23('zjwt',[3,0.0005],0)

若想看图中点的坐标可执行下面的命令：[t,y]= ode23('zjwt',[3,0.0005],0) plot(t,y)此时缉私艇的位置坐标是(0.00050000000000,1.96013657712118)执行下面的命令:

ode45('zjwt',[3,0.0005],0)

若想看图中点的坐标可执行下面的命令：[t,y]=ode45('zjwt',[3,0.0005],0) plot(t,y)此时缉私艇的位置坐标是(0.0005,1.9675 )

计算机仿真运算法

4计算机仿真法当建立动态系统的微分方程模型很困难时，我们可以用计算机仿真法对系统进行分析研究。所谓计算机仿真就是利用计算机对实际动态系统的结构和行为进行编程、模拟和计算，以此来预测系统的行为效果。

计算机仿真运算法

计算机仿真法缉私艇的初始位在点(c,0), (0, at k ) t t k:走私船的位置： ( x k, yk )缉私艇的位置：

y

走私船初始位在点(0,0),方向为y轴正方向，o(0, at k 1 )

c

x

追赶方向可用方向余弦表示为： (0, at k ) 0 xk k cos k (0 xk )2 (atk yk )2 ( xk, yk ) atk yk o sin k (0 xk )2 (atk yk )2时间步长

( xk 1, yk 1 ).

( xk 1, yk 1 ).则 t t k 1 t k t时,缉私艇的位置： xk 1 xk xk b t cos k, yk 1 yk yk b t sin k

计算机仿真运算法

仿真算法：第一步：设置时间步长 t,速度a, b及初始位置 x 0 c, y 0, k 0第二步：计算动点缉私艇D在时刻 t k 1 t k t时的坐标 xk xk 1 x k b t 2 2 xk (at k yk ) at yk y k 1 y k b t 2 xk (at k yk )2~~ t t t ( x, y )时的坐标计算走私船R在时刻k 1 k k 1 k 1

~ 0 x k 1~ yk 1 a ( t k t )

计算机仿真运算法

第三步：计算缉私艇与走私船这两个动点之间的距离：dk 2~ )2 ( y ~ ( x k 1 x y ) k 1 k 1 k 1

根据事先给定的距离，判断缉私艇是否已经追上了走私船，从而判断退出循环还是让时间产生一个步长，返回到第二步继续进入下一次循环；第四步：当从上述循环退出后，由点列 ( xk 1, yk 1 )~,~ ( x和 k 1 yk 1 )可分别绘制成两条曲线即为缉私艇和走私船走过的轨迹曲线。

计算机仿真运算法

显示船

与艇行进路线程序

c=3; a=0.4/60; b=0.8/60; d=0.01;dt=2;t=0; jstx=c;jsty=0;zscx=0;zscy=0; while (sqrt((jstx-zscx)^2+(jsty-zscy)^2)>d) pause(0.1) hold on axis([0,3,0,2]) t=t+dt; jstx=jstx-b*dt*jstx/sqrt(jstx^2+(a*t-jsty)^2); jsty=jsty+b*dt*(a*t-jsty)/sqrt(jstx^2+(a*tjsty)^2); zscy=a*t; plot(jstx,jsty,'rO',zscx,zscy, 'b*') end