试编程模拟导d打击敌机的动态过程，并实时给出飞机和导d的位置坐标。

%{

当t =0 时,导d位于原点（0，0）,敌机位于(0,100)点；

当时刻t ,导d位于L(x(t),y(t))，敌机位于(vt,100)点。 假设敌机的速度为v

导d速度可由水平分速度与垂直分速度合成：

(dx/dt)^2+(dy/dt)^2=(2v)^2______1

导d方向指向敌机,导d轨迹的导数就是其切线，所以

dy/dx=(100-y)/(vt-x)__________2

而dy/dx＝(dy/dt)/(dx/dt)

解以上微分方程组，初始条件为：x(0)=0,y(0)=0

用数值解法，可以用差分方程法，也可以用龙格库塔法，还可以消去t，化为二阶微分方程。这里用差分方法。

dx=x(k+1)-x(k);dy=y(k+1)-y(k);dt=t(k+1)-t(k)=h

%}

%Matlab程序：

clear;clc

h=01;%时间步长

k=1;

v=1;%假设敌机的速度为1，其它数值也行，不影响结果，只是步长h做相应的变化

t(1)=0;x(1)=0;y(1)=0;%初始值

while y<=100

x(k+1)=x(k)+2vh/sqrt(1+((100-y(k))/(vt(k)-x(k)))^2);

y(k+1)=y(k)+2vh/sqrt(1+((vt(k)-x(k))/(100-y(k)))^2);

t(k+1)=hk;

k=k+1;

end

plot(x,y,x(1):005:x(end),100)

t=t(end),x=x(end),y=y(end)

结果：

t = 667000

x = 666999

y = 1000000

导d总体优化软件是用于对导d进行设计和优化的软件，其功能包括对导d的结构、控制、导引、d头等方面进行优化。以下是一些常用的导d总体优化软件：

ABAQUS：ABAQUS是一个基于有限元分析（FEA）的软件平台，可用于模拟和优化导d的结构、d头和其他组件。

MATLAB：MATLAB是一种用于科学计算和工程应用的高级编程语言和交互式环境，可用于导d的动力学、控制和优化等方面。

MSC Nastran：MSC Nastran是一个FEA软件包，可用于导d的结构优化、疲劳分析和振动分析等方面。

Ansys：Ansys是一种FEA软件，可用于导d的结构、控制、传感器和d头等方面的优化。

Simulink：Simulink是MATLAB的一个扩展，可用于导d的动态建模、仿真和控制系统设计。

需要注意的是，这些软件都需要专业知识和技能才能正确使用。在使用这些软件之前，建议您先进行相关培训或获取相关领域的知识背景。

下面是camshift物体追踪代码，需要你用avi视频测试，matlab对avi视频格式要求比较严格。但是可以试试mmreader函数读取视频。

% Adam Kukucka

% Zach Clay

% Marcelo Molina

% CSE 486 Project 3

function [ trackmov probmov centers ] = camshift

%

% initialize vari ables

%

rmin = 0; %min row value for search window

rmax = 0; %max row value for search window

cmin = 0; %min col value for search window

cmax = 0; %max col value for search window

numofframes = 0; %number of frames in the avi

threshold = 1; %threshold for convergence

centerold = [0 0]; %for convergence previous center of window

centernew = [0 0]; %for convergence new center of window

%

% Pre code load movie and select initial frame

%

% prompt user for avi file name

%%%%%user_entry = input('Please enter an avi filename: ','s');

% load the avi file handle is M

%%%%M = aviread(user_entry);

M=aviread('8888avi');

% get number of frames

[dontneed numberofframes] = size(M)

% initialize matrix to hold center coordinates

imagecenters = zeros(numberofframes, 2);

% extract the first frame from the avi

Frame1 = M(1,1);

Image1 = frame2im(Frame1);

%%% images(:, :, numberofframes) = G(:,:);

% get search window for first frame

[ cmin, cmax, rmin, rmax ] = select( Image1 );

cmin = round(cmin);

cmax = round(cmax);

rmin = round(rmin);

rmax = round(rmax);

wsize(1) = abs(rmax - rmin);

wsize(2) = abs(cmax - cmin);

% create histogram

% translate to hsv

hsvimage = rgb2hsv(Image1);

% pull out the h

huenorm = hsvimage(:,:,1);

% scale to 0 to 255

hue = huenorm255;

% set unit type

hue=uint8(hue);

% Getting Histogram of Image:

histogram = zeros(256);

for i=rmin:rmax

for j=cmin:cmax

index = uint8(hue(i,j)+1);

%count number of each pixel

histogram(index) = histogram(index) + 1;

end

end

%

% Algorithm from pdf

%

aviobj1 = avifile('example3avi');

aviobj2 = avifile('example4avi');

% for each frame

for i = 1:200

disp('Processing frame');

disp(i);

Frame = M(1, i);

I = frame2im(Frame);

% translate to hsv

hsvimage = rgb2hsv(I);

% pull out the h

huenorm = hsvimage(:,:,1);

% scale to 0 to 255

hue = huenorm255;

% set unit type

hue=uint8(hue);

[rows cols] = size(hue);

% choose initial search window

% the search window is (cmin, rmin) to (cmax, rmax)

% create a probability map

probmap = zeros(rows, cols);

for r=1:rows

for c=1:cols

if(hue(r,c) ~= 0)

probmap(r,c)= histogram(hue(r,c));

end

end

end

probmap = probmap/max(max(probmap));

probmap = probmap255;

count = 0;

rowcenter = 0; % any number just so it runs through at least twice

colcenter = 0;

rowcenterold = 30;

colcenterold = 30;

% Mean Shift for 15 iterations or until convergence(the center doesnt

% change)

while (((abs(rowcenter - rowcenterold) > 2) && (abs(colcenter - colcenterold) > 2)) || (count < 15) )

%for j = 1:5

%disp('meanshift');

% disp(j);

rmin = rmin - 7; %increase window size and check for center

rmax = rmax + 7;

cmin = cmin - 7;

cmax = cmax + 7;

rowcenterold = rowcenter; %save old center for convergence check

colcenterold = colcenter;

[ rowcenter colcenter M00 ] = meanshift(I, rmin, rmax, cmin,

cmax, probmap);

% given image (I), search window(rmin rmax cmin cmax)

% returns new center (colcenter, rowcenter) for window and

% zeroth moment (Moo)

% redetermine window around new center

rmin = round(rowcenter - wsize(1)/2);

rmax = round(rowcenter + wsize(1)/2);

cmin = round(colcenter - wsize(2)/2);

cmax = round(colcenter + wsize(2)/2);

wsize(1) = abs(rmax - rmin);

wsize(2) = abs(cmax - cmin);

count = count + 1;

end

% mark center on image

%save image

G = 2989I(:,:,1)

+5870I(:,:,2)

+1140I(:,:,3);

trackim=G;

%make box of current search window on saved image

for r= rmin:rmax

trackim(r, cmin) = 255;

trackim(r, cmax) = 255;

end

for c= cmin:cmax

trackim(rmin, c) = 255;

trackim(rmax, c) = 255;

end

aviobj1 = addframe(aviobj1,trackim);

aviobj2 = addframe(aviobj2,probmap);

%create image movie, and probability map movie

trackmov(:,:,i)= trackim(:,:);

probmov(:,:,i) = probmap(:,:);

% save center coordinates as an x, y by doing col, row

centers(i,:) = [colcenter rowcenter];

% Set window size = 2 (Moo/256)^1/2

windowsize = 2 (M00/256)^5;

% get side length window size is an area so sqrt(Area)=sidelength

sidelength = sqrt(windowsize);

% determine rmin, rmax, cmin, cmax

rmin = round(rowcenter-sidelength/2);

rmax = round(rowcenter+sidelength/2);

cmin = round(colcenter-sidelength/2);

cmax = round(colcenter+sidelength/2);

wsize(1) = abs(rmax - rmin);

wsize(2) = abs(cmax - cmin);

end

% end for loop

% Adam Kukucka

% Zach Clay

% Marcelo Molina

% CSE 486 Project 3

function [ rowcenter colcenter M00 ] = meanshift(I, rmin, rmax, cmin,

cmax, probmap)

%inputs

% rmin, rmax, cmin, cmax are the coordiantes of the window

% I is the image

%outputs

% colcenter rowcenter are the new center coordinates

% Moo is the zeroth mean

%

% initialize

%

M00 = 0; %zeroth mean

M10 = 0; %first moment for x

M01 = 0; %first moment for y

histdim = (0:1:255); % dimensions of histogram 0 to 255, increment by 1

[rows cols] = size(I);

cols = cols/3; % 8

%

% Main code

%

% determine zeroth moment

for c = cmin:cmax

for r = rmin:rmax

M00 = M00 + probmap(r, c);

end

end

% determine first moment for x(col) and y(row)

for c = cmin:cmax

for r = rmin:rmax

M10 = M10 + cprobmap(r,c);

M01 = M01 + rprobmap(r,c);

end

end

% determine new centroid

% x is cols

colcenter = M10/M00;

% y is rows

rowcenter = M01/M00;

% Adam Kukucka

% Zach Clay

% Marcelo Molina

% CSE 486 Project 3

function [ cmin, cmax, rmin, rmax ] = select( I )

%UNTITLED1 Summary of this function goes here

% Detailed explanation goes here

% for array x is cols, y is rows

image(I);

k = waitforbuttonpress;

point1 = get(gca,'CurrentPoint'); %mouse pressed

rectregion = rbbox;

point2 = get(gca,'CurrentPoint');

point1 = point1(1,1:2); % extract col/row min and maxs

point2 = point2(1,1:2);

lowerleft = min(point1, point2);

upperright = max(point1, point2);

cmin = round(lowerleft(1));

cmax = round(upperright(1));

rmin = round(lowerleft(2));

rmax = round(upperright(2));

循环之前令sk=[];

在循环结构里，xk的值运算出来之后，

sk=[sk,xk]; %这样得出sk是行向量。

或者

sk=[sk;xk]; %这样得出的sk是列向量。

如果数据行数实在太多，超出command window所显示的行数，

那就

>>more on

这样就分页显示了。

这个错误别人不好判断，因为：

1、程序需要读文件，你没有提供；

2、你贴的出错信息也不完整，究竟是在这段代码的哪一行出的错？

从你使用plotm的两行代码看，应该都没什么大问题。

请把文件 mrwan_earthquake_catalogue_140228txt 上传到网盘，再帮你分析。

或者，你也可以试着自己解决问题：

在出错的位置设置断点，跟踪程序的运行，看究竟出错时变量的值是什么，以及为什么会是这个值。找到原因了，一般来说，也就知道该怎么解决了。

目录 一．设计目的·····································3二．设计思想·····································3三．设计步骤·····································3四．心得体会·····································10 一．设计目的1、学会用Matlab求简单微分方程的解析解 2、学会用Matlab求微分方程的数值解 二．设计思想 用Matlab软件求常微分方程的数值解，完成导d追踪问题，慢跑者与狗，地中海鲨鱼问题的求解。 三．设计步骤1 微分方程的解析解 求微分方程（组）的解析解命令:dsolve(‘方程1’, ‘方程2’,…‘方程n’, ‘初始条件’, ‘自变量’)记号: 在表达微分方程时，用字母D表示求微分，D2、D3等表示求高阶微分任何D后所跟的字母为因变量，自变量可以指定或由系统规则选定为确省 2 微分方程的数值解（一）常微分方程数值解的定义在生产和科研中所处理的微分方程往往很复杂且大多得不出一般解。而在实际上对初值问题，一般是要求得到解在若干个点上满足规定精确度的近似值，或者得到一个满足精确度要求的便于计算的表达式。因此，研究常微分方程的数值解法是十分必要的。

（二）建立数值解法的一些途径

1、用差商代替导数 若步长h较小，则有 故有公式 此即欧拉法2、使用数值积分对方程y’=f(x,y), 两边由xi到xi+1积分，并利用梯形公式，有 故有公式： 实际应用时，与欧拉公式结合使用： 此即改进的欧拉法。 3、使用泰勒公式以此方法为基础，有龙格-库塔法、线性多步法等方法 4、数值公式的精度 当一个数值公式的截断误差可表示为O（hk+1）时（k为正整数，h为步长），称它是一个k阶公式。k越大，则数值公式的精度越高。�6�1 欧拉法是一阶公式，改进的欧拉法是二阶公式。�6�1 龙格-库塔法有二阶公式和四阶公式。�6�1 线性多步法有四阶阿达姆斯外插公式和内插公式 （三）用Matlab软件求常微分方程的数值解 [t，x]=solver（’f’,ts,x0,options）注意:1、在解n个未知函数的方程组时，x0和x均为n维向量，m-文件中的待解方程组应以x的分量形式写成 2、使用Matlab软件求数值解时，高阶微分方程必须等价地变换成一阶微分方程组 （四）数学建模实例1 导d追踪问题设位于坐标原点的甲舰向位于x轴上点A(1, 0)处的乙舰发射导d，导d头始终对准乙舰如果乙舰以最大的速度v0(是常数)沿平行于y轴的直线行驶，导d的速度是5v0，求导d运行的曲线方程又乙舰行驶多远时，导d将它击中？ 解法一（解析法）假设导d在t时刻的位置为P(x(t), y(t))，乙舰位于 由于导d头始终对准乙舰，故此时直线PQ就是导d的轨迹曲线弧OP在点P处的切线，即有 即 (1)又根据题意，弧OP的长度为 的5倍，即 (2)由(1),(2)消去t整理得模型:

解即为导d的运行轨迹: 当 时 ，即当乙舰航行到点 处时被导d击中 被击中时间为: 若v0=1, 则在t=021处被击中 结果如图： 解法二 (建立参数方程求数值解) 设时刻t乙舰的坐标为(X(t)，Y(t))，导d的坐标为(x(t)，y(t)) 解导d运动轨迹的参数方程建立m-文件eq2m如下： 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); 取t0=0，tf=2，建立主程序chase2m如下： [t,y]=ode45('eq2',[0 2],[0 0]); Y=0:001:2; plot(1,Y,'-')， hold on plot(y(:,1),y(:,2),'') 结果见图1 导d大致在（1，02）处击中乙舰，与前面的结论一致在chase2m中，按二分法逐步修改tf，即分别取tf=1，05,025,…,直到tf=021时，得图2 结论：时刻t=021时，导d在（1，021）处击中乙舰。 2地中海鲨鱼问题 意大利生物学家Ancona曾致力于鱼类种群相互制约关系的研究，他从第一次世界大战期间,地中海各港口捕获的几种鱼类捕获量百分比的资料中，发现鲨鱼等的比例有明显增加（见下表），而供其捕食的食用鱼的百分比却明显下降显然战争使捕鱼量下降，食用鱼增加，鲨鱼等也随之增加，但为何鲨鱼的比例大幅增加呢？ 他无法解释这个现象，于是求助于著名的意大利数学家VVolterra，希望建立一个食饵—捕食系统的数学模型，定量地回答这个问题，建立m-文件shierm如下： function dx=shier(t,x) dx=zeros(2,1); dx(1)=x(1)(1-01x(2)); dx(2)=x(2)(-05+002x(1));其次，建立主程序sharkm如下： [t,x]=ode45('shier',[0 15],[25 2]); plot(t,x(:,1),'-',t,x(:,2),'') plot(x(:,1),x(:,2))求解结果 左图反映了x1（t）与x2（t）的关系。 可以猜测： x1（t）与x2（t）都是周期函数。模型（二） 考虑人工捕获模型求解:1、分别用m-文件shier1m和shier2m定义上述两个方程2、建立主程序shark1m, 求解两个方程，并画出两种情况下鲨鱼数在鱼类总数中所占比例 x2(t)/[x1(t)+x2(t)] 图中实线为战前的鲨鱼比例，“”线为战争中的鲨鱼比例 结论：战争中鲨鱼的比例比战前高 （五）源程序代码 clearx=0:001:1;y=-5(1-x)^(4/5)/8+5(1-x)^(6/5)/12+5/24;plot(x,y,'') t0=0;tf=2;[t,y]=ode45('eq2',[t0 tf],[0 0]);X=1;Y=0:001:2;plot(X,Y,'-')hold onplot(y(:,1),y(:,2),'') t0=0;tf=375;[t,y]=ode45('eq3',[t0 tf],[0 0]);T=0:01:2pi;X=10+20cos(T);Y=20+15sin(T);plot(X,Y,'-')hold onplot(y(:,1),y(:,2),'') function dy=eq1(x,y) dy=zeros(2,1); dy(1)=y(2);dy(2)=1/5sqrt(1+y(1)^2)/(1-x); 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); dsolve('Du=1+u^2','t') y=dsolve('D2y+4Dy+29y=0','y(0)=0,Dy(0)=15','x')[x,y,z]=dsolve('Dx=2x-3y+3z','Dy=4x-5y+3z','Dz=4x-4y+2z','t'); x=simple(x) y=simple(y) z=simple(z) [T,Y]=ode15s('vdp1000',[0 3000],[2 0]); plot(T,Y(:,1),'-')[T,Y]=ode45('rigid',[0 12],[0 1 1]); plot(T,Y(:,1),'-',T,Y(:,2),'',T,Y(:,3),'+') x0=0;xf=099999;[x,y]=ode23('eq1',[x0 xf],[0 0]); Y=0:001:2; plot(1,Y,'g') hold onplot(x,y(:,1),'b') function dy=rigid(t,y) dy=zeros(3,1); dy(1)=y(2)y(3); dy(2)=-y(1)y(3); dy(3)=-051y(1)y(2); [t,x]=ode45('shier',[0 15],[25 2]); plot(t,x(:,1),'-',t,x(:,2),'') figure(2) plot(x(:,1),x(:,2)) [t1,x]=ode45('shier1',[0 15],[25 2]); [t2,y]=ode45('shier2',[0 15],[25 2]); x1=x(:,1);x2=x(:,2); x3=x2/(x1+x2); y1=y(:,1);y2=y(:,2); y3=y2/(y1+y2); plot(t1,x3,'-',t2,y3,'')function dx=shier(t,x) dx=zeros(2,1); dx(1)=x(1)(1-01x(2));dx(2)=x(2)(-05+002x(1)); function dx=shier1(t,x) dx=zeros(2,1); dx(1)=x(1)(07-01x(2));dx(2)=x(2)(-08+002x(1)); function dy=shier2(t,y) dy=zeros(2,1); dy(1)=y(1)(09-01y(2)); dy(2)=y(2)(-06+002y(1));function dy=vdp1000(t,y) dy=zeros(2,1); dy(1)=y(2); dy(2)=1000(1-y(1)^2)y(2)-y(1);四．心得体会通过本次Matlab的课程设计，我对Matlab基础知识有了深刻了解，了解了Matlab软件在数学建模中的使用。学会了用Matlab解简单微分方程的数值解和解析解。通过本次课程设计我还认识到了坚持和认真的重要性，虽然在任务的完成过程中遇到了很多的难题，但通过努力查阅书本和网络资料都一一解决了。本次课程设计是我受益良多，为我以后对Matlab的熟练使用打下了良好的基础。非常感谢学校安排这次课程设计。

以上就是关于试编程模拟导d打击敌机的动态过程，并实时给出飞机和导d的位置坐标。全部的内容，包括:试编程模拟导d打击敌机的动态过程，并实时给出飞机和导d的位置坐标。、导d总体优化软件有哪些、关于MATLAB中多目标的跟踪。等相关内容解答，如果想了解更多相关内容，可以关注我们，你们的支持是我们更新的动力！