代朋友请教如何提高程序的运行效率

20wangz

我的一个朋友，最近在做最优控制方面的学习和研究。他写了一个程序，但是运行效率特别低，运行时间经常以小时计。我看了一下他的程序，别的咱不懂，但程序中用了很多的循环，难怪快不起来。我向他介绍了仿真论坛上有许多热心的高手，于是他委托我把问题贴上来，向各位兄台请教：怎样才能让程序尽快地得到结果呢？

下面是他的程序代码：

1. clc
   
2. clear
   
3. %已知参数
   
4. h=-36.884;    %热机械参数1
   
5. s=-18.378;    %热机械参数2
   
6. a0=2*10^6;     %热机械参数3
   
7. K0=1;          %热机械参数4
   
8. T0=300;        %环境温度
   
9. p0=200;        %外部压力
   
10. R=8.3145;      %气体常数
    
11. Cv=2.5*R;      %热熔与R的关系
    
12. p1=16*p0/1000;      %无量纲系数p'
    
13. b1=0.4;           %牛顿传热规律时的无量纲系统的系数
    
14. b2=1;             %牛顿传热规律时的无量纲系统的系数
    
15. b3=0.52632;         %牛顿传热规律时的无量纲系统的系数
    
16. b4=98468;           %牛顿传热规律时的无量纲系统的系数
    
17. b5=3.8643*10^(-8);       %牛顿传热规律时的无量纲系统的系数
    
18. b6=1.0211;                 %牛顿传热规律时的无量纲系统的系数
    
19. b7=128.368;                  %牛顿传热规律时的无量纲系统的系数
    
20. t=0.5;                      %膨胀和压缩冲程的时间
    
21. x(1)=1;                   %无量纲位移初值
    
22. xfin=2;                    %无量纲位移末值
    
23. T(1)=1.6280;                %无量纲温度的初值
    
24. Tfin=1.2464;                   %无量纲温度的末值
    
25. num=200;
    
26. detat=t/(num-1);
    
27. ww(1)=0;
    
28. k=1;
    
29. for j=0.000001:0.001:1
    
30.      n=1;
    
31.      lamda1(1)=j;
    
32.      for jj=0.00001:0.001:6
    
33.          lamda2(1)=jj;
    
34.          v(1)=-1/(2*b6*x(1))*(b1*lamda1(1)*T(1)-lamda2(1)*x(1)-T(1));
    
35.          sigma(1)=exp((-1*h+s*T(1))/T(1));
    
36.          L(1)=(x(1)^2+b4*b5*x(1)*T(1)*sigma(1))^(1/2);
    
37.           for i=2:num
    
38.               ww(i)=(i-1)*detat;
    
39.               lamda2(i)=lamda2(i-1)+detat*(v(i-1)*T(i-1)/(x(i-1)^2)-lamda1(i-1)*(-b2*(T(i-1)^4-1)+b1*v(i-1)*T(i-1)/(x(i-1)^2)+a0*b3/(b4*T(i-1)*sigma(i-1))*exp(-1*a0/(b4*sigma(i-1)*T(i-1))*(-1*x(i-1)+L(i-1)))*(-1+(2*x(i-1)+b4*b5*sigma(i-1)*T(i-1))/(2*L(i-1)))));
    
40.               lamda1(i)=lamda1(i-1)+detat*(-1*v(i-1)/x(i-1)-lamda1(i-1)*(-1*b1*v(i-1)/x(i-1)-4*(1+b2*x(i-1))*T(i-1)^3-b3*exp(-1*a0/(b4*sigma(i-1)*T(i-1))*(-1*x(i-1)+L(i-1)))*a0/(b4*T(i-1)^2*sigma(i-1))*(-1*x(i-1)+L(i-1)+h*(-1*x(i-1)+L(i-1))/T(i-1)-b4*b5*x(i-1)*sigma(i-1)*(T(i-1)+h)/(2*L(i-1)))));
    
41.               T(i)=T(i-1)+detat*(-1*b1*T(i-1)*v(i-1)/x(i-1)+b3*(1-exp(-1*a0/(b4*T(i-1)*sigma(i-1))*(-1*x(i-1)+L(i-1))))-(1+b2*x(i-1))*(T(i-1)^4-1));
    
42.               x(i)=x(i-1)+detat*v(i-1);
    
43.               v(i)=-1/(2*b6*x(i))*(b1*lamda1(i)*T(i)-lamda2(i)*x(i)-T(i));
    
44.               sigma(i)=exp((-1*h+s*T(i))/T(i));
    
45.               L(i)=(x(i)^2+b4*b5*x(i)*T(i)*sigma(i))^(1/2);
    
46.           end
    
47.           zz(k,n)=T(num);
    
48.           zzz(k,n)=x(num);
    
49.           n=n+1;
    
50.      end
    
51.      k=k+1;
    
52. end
    
53.      
    
54. for i=1:k-1
    
55.     for j=1:n-1
    
56.         if ((abs(zz(i,j)-Tfin)<0.001)&&(abs(zzz(i,j)-xfin)<0.001))
    
57.             i
    
58.             j
    
59.         end
    
60.     end
    
61. end   
    
62.          

复制代码

有关的说明见附件。

代他先谢谢大家了！

pingnator

1# 20wangz


根据你的code 我把前面的2层for 循环用向量运算代替了。试了一下，按照你的参数设置，执行完大概26分钟左右。

但是有几个问题：
（1） 你的ww 赋值了，但是从来都没用；
（2） 你自己程序里面的44 行在运算过程中会出现singularity的情况， 因为T可能非常非常接近0; 应该根据这行代码的实际物理意义加以判断排除；


clc
clear
close

h=-36.884;   
s=-18.378;   
a0=2*10^6;     
K0=1;         
T0=300;        
p0=200;        
R=8.3145;      
Cv=2.5*R;      
p1=16*p0/1000;      
b1=0.4;           
b2=1;            
b3=0.52632;         
b4=98468;           
b5=3.8643*10^(-8);      
b6=1.0211;                 
b7=128.368;                  
t=0.5;                     
x0 =1;                     
xfin=2;                     
T0=1.6280;               
Tfin=1.2464;               
num=200;
detat=t/(num-1);
ww(1)=0;
k=1;

[lamda1,lamda2] = meshgrid(0.000001:0.001:1,0.00001:0.001:6);
x = x0* ones(size(lamda1));
T = T0* ones(size(lamda1));
v = -1./(2*b6*x).*(b1*lamda1.*T-lamda2.*x-T);
sigma=exp((-1*h+s*T)./T);
L=(x.^2+ b4*b5*x.*T.*sigma).^(1/2);
for i=2:num
    i
    lamda1_old = lamda1; lamda2_old = lamda2; T_old = T; x_old = x; v_old = v; sigma_old = sigma; L_old = L;
    lamda2=lamda2_old+detat*(v_old.*T_old./(x_old.^2)-lamda1_old.*(-b2*(T_old.^4-1)+b1*v_old.*T_old./(x_old.^2)+a0*b3./(b4*T_old.*sigma_old).*exp(-1*a0./(b4*sigma_old.*T_old).*(-1*x_old+L_old)).*(-1+(2*x_old+b4*b5*sigma_old.*T_old)./(2*L_old))));
    lamda1=lamda1_old+detat*(-1*v_old./x_old-lamda1_old.*(-1*b1*v_old./x_old-4*(1+b2*x_old).*T_old.^3-b3*exp(-1*a0./(b4*sigma_old.*T_old).*(-1*x_old+L_old))*a0./(b4*T_old.^2.*sigma_old).*(-1*x_old+L_old+h*(-1*x_old+L_old)./T_old-b4*b5*x_old.*sigma_old.*(T_old+h)./(2*L_old))));
    T = T_old+detat*(-1*b1*T_old.*v_old./x_old+b3*(1-exp(-1*a0./(b4*T_old.*sigma_old).*(-1*x_old+L_old)))-(1+b2*x_old).*(T_old.^4-1));
    x=x_old+detat*v_old;
    v=-1./(2*b6*x).*(b1*lamda1.*T-lamda2.*x-T);
    sigma = exp((-1*h+s*T)./T);
    L =(x.^2+b4*b5*x.*T.*sigma).^(1/2);
end
zz = T;
zzz = x;
[Nrow, Ncol]= find((abs(zz-Tfin)<0.001)&&(abs(zzz-xfin)<0.001));

20wangz

刚才给朋友打了电话，他非常高兴！让我替他再次感谢回答问题的兄台。

谢谢！