matlab 内弹道计算问题

　　MATLAB是一种对技术计算高性能的语言。
　　它集成了计算，可视化和编程于一个易用的环境中，在此环境下，问题和解答都表达为我们熟悉的数学符号。
　　典型的应用有：

数学和计算
算法开发
建模，模拟和原形化
数据分析，探索和可视化
科学与工程制图
应用开发，包括图形用户界面的建立
MATLAB是一个交互式的系统，其基本数据元素是无须定义维数的数组。
　　这让你能解决很多技术计算的问题，尤其是那些要用到矩阵和向量表达式的问题。而要花的时间则只是用一种标量非交互语言(例如C或Fortran)写一个程序的时间的一小部分。

一、字符和字符串
在MATLAB中不管是字符，还是字符串都是用单引号。而在C、C++、C#、Java等语言中，字符是用单引号的，字符串则必须用双引号。因此在MATLAB中如果需要在字符串中用到单引号的话，就要用两个单引号来代替。
如：
>> a='this''is an apple'
a =
this'is an apple
同理，如果字符串中需要双引号的话，可以用四个单引号来代替。如下：
>> b='I''''am Tim'
b =
I''am Tim
在MATLAB中的字符串连接，我们可以使用一对中括号。当然，这种连接方式也适用于向量、矩阵等的连接。如下：
>> c=[a b]
c =
this'is an appleI''am Tim
>> a1=[1 2];b1=[3 4];
>> c1=[a1 b1]
c1 =
1 2 3 4
>> c2=[a1;b1]
c2 =
1 2
3 4
但是，对于字符串的连接不能使用加号（+）来进行。这点和在C++、C#、Java等语言中是不一样的。因为在MATLAB中这些字符串也是以矩阵的形式存储的，你要是用加号的话，需要两个矩阵的大小一致。比如：
>> a2='hello';b2='mustb';
>> a2+b2
ans =
213 218 223 224 209
但是，很明显，加完之后都是一系列的值了。本来，这些字符串也是以数值的形式存储的。而我们要取的字符串中某一个字符的时候，也是很方便的，直接像引用矩阵的元素一样。如下：
>> a2(1)
ans =
h
>> a2(4)
ans =
l
>> find(a2>'i')
ans =
3 4 5
可以看到，a2里面在字母i后面的字符有第3、4、5个字符，也就是llo了。
在MATLAB中要将一个数字转换为字符可以用num2str，将字符转换为数字可以用str2num。如下：
>> stra1=num2str(a1)
stra1 =
1 2
>> class(stra1)
ans =
char
>> strb1='[3 4]';
>> str2num(strb1)
ans =
3 4
对于strb1，我们可以使用eval函数，该函数可以将字符串作为一个MATLAB命令去执行。如下：
>> p=eval(strb1)
p =
3 4
相当于p=[3 4]。
当然了，C语言里面的sprintf、fprintf函数我们都是可以在MATLAB中使用的。比如：
>> sprintf('Tim is %d years old.',24)
ans =
Tim is 24 years old.
至于fprintf函数，则需要用fopen去打开一个文件，然后写入数据进去。如下：
>> fd=fopen('a.txt','w+');
>> fprintf(fd,'I am Tim.
');
>> fprintf(fd,'My age is %d
',24);
>> fclose(fd);
然后，用MATLAB Editor打开MATLAB当前目录下的a.txt就可以看到：

但是，当你用记事本打开的时候会看到：

这是因为在Windows上要想换行需要
，而不能单是
。

二、if语句、for语句、switch语句、while语句、try-catch语句都以end结尾
在MATLAB中，没有goto和do…loop和do…while语句。而且语句的结尾都是以end结尾的。
>> t=4;
>> if(t<5),disp('t is smaller than 5.'),end
t is smaller than 5.
当然，if语句也可以有else语句。
>> if(t>5),disp('t is bigger than 5.'),else,disp('t is smaller than 5.'),end
t is smaller than 5.
还可以用elseif。
>> if(t>5),disp('t is bigger than 5.'),elseif(t==5),disp('t is equal 5.'),else,disp('t is smaller than 5.'),end
t is smaller than 5.
在MATLAB中，由于矩阵的索引都是从1开始的，这点和C、C++、C#、JAVA等语言是不一样的。所以，在循环的时候建议从1开始循环，而且循环的时候不推荐用变量i和j。因为这两个变量在MATLAB中表示的是复数的虚部变量。当然，如果你的程序里面没有复数，就不用怕。如果有复数，就要小心了。
之所以从1开始，是因为大多数时候我们使用循环是为了依次能访问到矩阵的每一个元素，所以，如果矩阵的索引出现了0，那么就会报出错误。所以，从1开始循环，是一种很好的习惯，而不是要求你必须这样做。
而且，在for循环中，你的循环变量不仅可以是一个数，也可以是一个向量。看下面的程序：
>> aa=[1 2;3 4; 5 6;7 8];
>> for i=aa,disp(i),end
1
3
5
7
2
4
6
8
>> aa
aa =
1 2
3 4
5 6
7 8
也就是说，i第一次的值是aa的第一列，i第二次的值是aa的第二列。这样的话，这个循环变量i就不仅仅是一个数了，二是一列数。
在MATLAB中，switch语句中是不用break的，而且默认的选项不是default，而是otherwise。如下：
a=5;
switch(a),
case 1,
disp('a is 1.');
case 2,
disp('a is 2.');
case 5,
disp('a is 5.');
otherwise,
disp('sorry, i do not know.');
end
还需要注意的是case后面没有分号，而C、C++、C#、Java等语言中是必须是分号的。在MATLAB中可是没有分号，我程序中的逗号，也可以是没有的。指向上面的程序，在MATLAB中的输出是：
a is 5.
至于while语句的使用和C、C++、C#、Java等语言中的while基本上是一致的。如下：
a=5;
b=0;
while a>0
b=b+a;
a=a-1;
end
b
a
Matlab命令窗口的输出是这样的：
b =
15
a =
0
try-catch语句一般是用来捕捉错误的。
try
a=input('Input a number: ');
catch
error('unknown error.');
end
当运行上面的程序，在MATLAB中输入：
Input a number: #
??? Error using ==> tim_try at 4
unknown error.

三、变量不用声明
在MATLAB中的变量是可以不用声明的，这点和C类语言和Java是有很大区别的。但是，为了程序更清晰和加快程序的执行速度，提供初始化和声明有时是必要的。不如，
b=[];
for i=1:10
b=[b i];
end
b
在MATLAB命令窗口的输出如下：
b =
1 2 3 4 5 6 7 8 9 10
如果你不初始化b为空的话，如下：
for i=1:10
b=[b i];
end
b
在命令窗口的输出如下：??? Undefined function or variable 'b'.
上面说过了，中括号[]可以用来合并向量的，这里变量b不存在，因此不能合并。所以会出错。

四、没有++、--操作符
在MATLAB中是没有++和--操作符的，因此需要转换。如a++; 可以转化为a=a+1;至于++a；你只需要将a=a+1的位置放到前面就可以了。

五、点乘、点除，但没有点加、点减
在点乘，是矩阵中的对应元素相乘。而不是矩阵的乘法。点除也是同样道理，对应元素想除。至于点加，就相当于加法。因为加减法本身就是对应元素的加减。

六、左除和右除
在C类和Java等语言中，除法就是一个操作符(/)，但在MATLAB中(/)和(\)是有区别的。/代表的是右除，\代表的是左除。
>> A=[1 2;3 4];B=[5 6;7 8];
>> C1=A\B
C1 =
-3 -4
4 5
所以A左除B，就相当于A的逆矩阵乘B。
>> inv(A)*B
ans =
-3.0000 -4.0000
4.0000 5.0000
那么可想而知，A右除B，就相当于A乘B的逆矩阵。
>> A/B
ans =
3.0000 -2.0000
2.0000 -1.0000
>> A*inv(B)
ans =
3.0000 -2.0000
2.0000 -1.0000

七、MATLAB的编程思想
当你用MATLAB来写程序的时候，尽量以向量、矩阵为单位来考虑问题。也就是在MATLAB中常常能用一个点乘来代替多个循环的操作。也就是用向量化的操作来代替循环，这样能够大大的提高MATLAB程序的运行速度。
举个简单的例子：
clear all
clc
a=round(rand(100)*100);
b=round(rand(100)*100);

tic
f1=a.*b;
toc

tic
for i=1:100,
for j=1:100,
f2(i,j)=a(i,j)*b(i,j);
end
end
toc
在命令窗口的输出如下：
Elapsed time is 0.000035 seconds.
Elapsed time is 0.001871 seconds.
如果我们给f2初始化下，看看能不能加快速度。
clear all
clc
a=round(rand(100)*100);
b=round(rand(100)*100);

tic
f1=a.*b;
toc

f2=zeros(100);
tic
for i=1:100,
for j=1:100,
f2(i,j)=a(i,j)*b(i,j);
end
end
toc
此时，命令窗口的输出如下：
Elapsed time is 0.000024 seconds.
Elapsed time is 0.000147 seconds.
可以看到后面的循环的时间明显提高了很多。那么我们为上面的点乘运算的f1也初始化的话，速度会如何？
Elapsed time is 0.000022 seconds.
Elapsed time is 0.000146 seconds.
可以看到速度提升的不是很明显。如果你将a和b的维数增长到10000的话，提示的速度就明显了很多了。

原因

所给链接的代码不完整，缺少ndd_fun函数。

顺便鄙视一下该链接的提供者，这么广为流传的东西下载居然还要积分，简直穷疯了。

代码

帮你好好找了一下，找到了完整的程序，供参考（全部代码保存到一个M文件运行即可，或直接下载附件）：

function ndd
%59nian130
A=0.87;         %枪(炮)膛横断面积A  dm^2
G=19;%33.4;           %弹重  kg
W0=2.04;        %药室容积  dm^3
l_g=25.0;       %身管行程  dm
P_0 =30000;       %起动压力  kpa
fai1=1.02;        %次要功系数
K=1.03;        %运动阻力系数φ1
theta =0.2;       %火药热力系数
%=========================================
f=950000;          %火药力  kg*dm/kg 
alpha=1;         %余容  dm^3/kg
delta=1.6;         %火药重度γ
%==================================
ome=2.2;%12.9;      %第一种装药量  kg
u1=5.0024*10^-5;        %第一种装药烧速系数  dm^3/(s*kg)
n1=0.82;         %第一种装药的压力指数n1
lambda=-0.0071;     %第一种装药形状特征量λ1
lambda_s=0;     %第一种装药分裂点形状特征量λ1s
chi=1.00716;          %第一种装药形状特征量χ1
chi_s=0;          %第一种装药分裂点形状特征量χ1s
mu=0;            %第一种装药形状特征量μ1
et1=1.14*10^-2;            %第一种装药药厚δ01
d1=2.5*10^-2;             %第一种装药火药内径d1
Ro1=0;             %药型系数α1
%=========================================
%常数与初值计算-----------------------------------------------------
l_0=W0/A;
Delta=ome/W0;
phi=K + ome/(3*G);
v_j=196*f*ome/(phi*theta*G);
v_j=sqrt(v_j);
B = 98*(et1*A)^2/( u1*u1*f*ome*phi*G );
B=B*(f*Delta)^(2-2*n1);
Z_s=1+Ro1*(d1/2+et1)/et1;
p_0=P_0/(f*Delta);
psi_0=(1/Delta - 1/delta)/(f/P_0 + alpha - 1/delta);
Z_0=(sqrt(1+4*psi_0*lambda/chi) - 1)/(2*lambda);
%解算子------------------------------------------------------------
C = zeros(1,12);
C(1)=chi;C(2)=lambda;C(3)=lambda_s;C(4)=chi_s;C(5)=Z_s;%
C(6)=theta;C(7)=B;C(8)=n1;C(9)=Delta;C(10)=delta;C(11)=alpha;C(12)=mu;
C;
y0=[Z_0;0;0;psi_0];
options = odeset('outputfcn','odeplot');
[tt,y] = ode45(@ndd_fun,0:100,[Z_0;0;0],options,C);
l = y(:,2);
l = l*l_0;
fl = find(l>=l_g);
fl = min(fl);
[tt,y] = ode45(@ndd_fun,0:0.005:fl,[Z_0;0;0],options,C);
Z = y(:,1);lx = y(:,2); vx = y(:,3); 
psi = (Z>=0&Z<1).*( chi*Z.*(1 + lambda*Z + mu*Z) ) +...%%%%%%%%%
      (Z>=1&Z<Z_s).*( chi_s*Z.*(1 + lambda_s*Z) ) +...
      (Z>=Z_s)*1;
l_psi = 1 - (Delta/delta)*(1-psi) - alpha*Delta*psi;
px = ( psi - vx.*vx )./( lx + l_psi );
p = px*f*Delta/100;
v = vx*v_j/10;
l = lx*l_0;
t = tt*l_0*1000/v_j;
fl = find(l>=l_g);
fl = min(fl)+1;
p(fl:end)=[];v(fl:end)=[];l(fl:end)=[];t(fl:end)=[];
pd=px*f*Delta/100/(1+ome/3/fai1/G);
pt=pd*(1+ome/2/fai1/G);
aa=max(px);
M=find(px==aa);
Pm=[tt(M)*l_0*1000/v_j lx(M)*l_0 vx(M)*v_j/10 px(M)*f*Delta/100 pt(M) pd(M) psi(M) Z(M)];
%ll=length(tt);
ran=find(Z>=1);
ran=min(ran);
Zf=[tt(ran)*l_0*1000/v_j lx(ran)*l_0 vx(ran)*v_j/10 px(ran)*f*Delta/100 pt(ran) pd(ran) psi(ran) Z(ran)];
jie=find(psi>=1);
jie=min(jie);
psij=[tt(jie)*l_0*1000/v_j lx(jie)*l_0 vx(jie)*v_j/10 px(jie)*f*Delta/100 pt(jie) pd(jie) psi(jie) Z(jie)];
pg=[tt(end)*l_0*1000/v_j lx(end)*l_0 vx(end)*v_j/10 px(end)*f*Delta/100 pt(end) pd(end) psi(end) Z(end)];
Ry1=[Zf;psij;pg;Pm];
Ry2=[tt*l_0*1000/v_j lx*l_0 vx*v_j/10 px*f*Delta/100 pt pd psi Z];
subplot(2,2,1);
plot(t,p,'linewidth',2);
grid on;
xlabel('\fontsize{8}\bft  (ms)');
ylabel('\fontsize{8}\bfp  (kg/cm^{2})');
title('\fontsize{8}\bft-p曲线');
subplot(2,2,2)
plot(t,v,'linewidth',2);
grid on;
xlabel('\fontsize{8}\bft  (ms)');
ylabel('\fontsize{8}\bfv  (m/s)');
title('\fontsize{8}\bft-v曲线');
subplot(2,2,3)
plot(l,p,'linewidth',2);
grid on;
xlabel('\fontsize{8}\bfl  (dm)');
ylabel('\fontsize{8}\bfp  (kg/cm^{2})');
title('\fontsize{8}\bfl-p曲线');
subplot(2,2,4)
plot(l,v,'linewidth',2);
grid on;
xlabel('\fontsize{8}\bfl  (dm)');
ylabel('\fontsize{8}\bfv  (m/s)');
title('\fontsize{8}\bfl-v曲线');
tspan = length(t)/20;
tspan = 1:ceil(tspan):length(t);
tspan(end) = length(t);
fprintf('        t(ms)     p(kg/cm^2)     v(m/s)       l(dm)');
format short g;
Result = [t(tspan) p(tspan) v(tspan) l(tspan)]
format;

%------------------------------------------------------------------

function dy = ndd_fun(t,y,C)
chi=C(1);lambda=C(2);lambda_s=C(3);chi_s=C(4);Z_s=C(5);mu=C(12);
theta=C(6);B=C(7);V=C(8);Delta=C(9);delta=C(10);alpha=C(11);
Z = y(1); l = y(2); v = y(3);
psi = (Z>=0&Z<1).*( chi*Z.*(1 + lambda*Z + mu*Z) ) +...
(Z>=1&Z<Z_s).*( chi_s*Z.*(1 + lambda_s*Z) ) +...
(Z>=Z_s)*1;
l_psi = 1 - (Delta/delta)*(1-psi) - alpha*Delta*psi;
p = ( psi - v*v )/( l + l_psi );
dy(1) = sqrt(theta/(2*B))*(p^V)*(Z>=0&Z<=Z_s);
dy(2) = v;
dy(3) = theta*p/2;
dy = [dy(1);dy(2);dy(3)];

结果