How to vectorize this code?

1 次查看（过去 30 天）

Taha Abbas Bin Rashid 2019-3-13

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/449740-how-to-vectorize-this-code

编辑： per isakson 2019-3-13

clear
clc
tic
%concentration
%[M]=x1 [C]=x2 [Dr]=x3 [Cx]=x4 [R^o]=x5 [Rm^o]=x6 [Pr]=x7
x=[100 1 1 0 0 0 0];
v=1e-18;  % control volume
NA=6.022e23; % Avogadro's Number
t=0;
nmax=2000;
%number of molecules
n=(NA*v).*x;
N=n;
%kp=k1  ktc0=k2  ktd0=k3  ktr=k4   ka=k5 kd=k5
k=[1516 3.469e8 0 0.22 0.45 1.1e7];
kmc=k./((NA*v));
kmc2=(2.*k)./((NA*v));
iter=1;
iter_t=0;
for n=1:nmax
    for t=t
        R(1)=kmc2(5)*N(3)*N(2);
        R(2)=kmc2(6)*N(5)*N(4);
        R(3)=kmc2(1)*N(5)*N(1);
        R(4)=kmc2(4)*N(5)*N(1);
        R(5)=kmc2(3)*N(5)*N(6)+kmc(4)*N(5)*N(6);
        Rt=sum(R);
        P=[R(1)/Rt  (R(1)+R(2))/Rt (R(1)+R(2)+R(3))/Rt (R(1)+R(2)+R(3)+R(4))/Rt (R(1)+R(2)+R(3)+R(4)+R(5))/Rt];
        Pt=sum(P);
        r1=rand;
        r2=rand;
        tau=(1/Rt)*log(1/r2);
        meu=[1:1:5];
        if  r1<Pt
            %r3=rand;
            if (0<r1)&(r1<=P(1))
                
                N(2)=N(2)-1;
                N(3)=N(3)-1;
                N(4)=N(4)+1;
                N(5)=N(5)+1;
                
                
                % r4=rand
            elseif (P(1)<r1)&(r1<=P(2))
                
                N(5)=N(5)-1;
                N(2)=N(2)+1;
                N(4)=N(4)-1;
                N(3)=N(3)+1;
                
                % r5=rand
            elseif (P(2)<r1)&(r1<=P(3))
                
                N(1)=N(1)-1;
                N(5)=N(5)-1;
                N(6)=N(6)+1;
                
                % r6=rand
            elseif (P(3)<r1)&(r1<=P(4))
                
                N(1)=N(1)-1;
                N(7)=N(7)+1;
                
                %r7=rand
            elseif (P(4)<r1)&(r1<=P(5))
                
                N(5)=N(5)-2;
                N(7)=N(7)+2;
            end
            NN(iter,:)=N;
            N(1,:)=NN(iter,:);
            %
            
        else
            "mistake"
        end
        % iter=iter+1;
        %NN(iter,:)=N;
        %plot(t_save,NN(:,2))
    end
    %y=N(2)./(NA*v)
    NN(n,:)=N;
    N(1,:)=NN(n,:);
    Rt_save(n,:)=Rt;
    
    t_save(n,:) =t;
    tau_save(n,:)=tau;
    y=(N(1)./NN(:,1));
    y=vpa(y);
    z=100*(NN(:,1)/N(1));
    figure (1)
    plot(t_save,y,'-s')
    xlabel('polymerization Time(s)');
    ylabel('(M0/M)')
    figure (2)
    plot(t_save,z,'-o')
    xlabel('polymerization Time(s)');
    ylabel('Conversion')
    % hold on
    t=t+tau;
    
end
toc

4 个评论
显示 2更早的评论隐藏 2更早的评论

Geoff Hayes 2019-3-13

Try preallocating or pre-sizing those arrays that you update on each iteration of the loop (this would be NN and N for example). See preallocation for details.

Taha Abbas Bin Rashid 2019-3-13