How can I speed up my Code?

Question

Gerrit 2020-1-9

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/499639-how-can-i-speed-up-my-code

编辑： Gerrit 2020-1-11

在 MATLAB Online 中打开

Good evening,

I have a problem with the performance of my Code, as Im new to matlab.

I hope the comments help you to understand everything.

The main problem is, that most of the time I have no points at the x-coordinate. Thats why I interpolate the closest and the secound closest. But that I have to do for every AeroDataId, so every Mach/CL combination in my database. To speed it up a little bit, I only use the points in a certin mach and Cl delta.

load('new_upper_airfoil');
%% searched point
Ma=0.7;
C_l=0.65;
%% delta of cl and Mach
cl_delta=0.015;
Ma_delta=0.015;
%% delta of the points in the original pressure distribution
Delta=2/height(originalDV);
DV=zeros(height(originalDV),2);
Temp=array2table(zeros(0,4),'VariableNames',{'X_C','CP','Ma','CL_res'});
Temp_update=array2table(zeros(0,4),'VariableNames',{'X_C','CP','Ma','CL_res'});
i=1;
%% load only the points within a Mach Number and Lift coefficient delta
new_upper_airfoil=new_upper_airfoil(new_upper_airfoil.Ma < (Ma+Ma_delta) & new_upper_airfoil.Ma > (Ma-Ma_delta) & new_upper_airfoil.CL_res < (C_l+cl_delta) & new_upper_airfoil.CL_res > (C_l-cl_delta),:);
new_lower_airfoil=new_lower_airfoil(new_lower_airfoil.Ma < (Ma+Ma_delta) & new_lower_airfoil.Ma > (Ma-Ma_delta) & new_lower_airfoil.CL_res < (C_l+cl_delta) & new_lower_airfoil.CL_res > (C_l-cl_delta),:);
%% Interpolation of upper_airfoil
for x_c=0:Delta:1
    %% Vector with all unique Ids (muss nicht nochmal für Unterseite erstellt werden)
    DataId_unique=unique(new_upper_airfoil.AeroDataId); %könnte gespeichert werden für alle Daten
   
    for a=1:1:size(DataId_unique)
        %% temp table with all the information for one ID (AeroData ID is no longer needed)
        new_upper_airfoil_unique=new_upper_airfoil(new_upper_airfoil.AeroDataId==DataId_unique(a,:),2:5); % könnte später gespeichert werden
        
        %% try to find x/c
        FIND=new_upper_airfoil_unique(new_upper_airfoil_unique.X_C==x_c,:);    
        
        if(isempty(FIND))
            %% find the closest and secound closest x-coordinate and interpolate between them
            [d,sortIndex]=sort(abs(new_upper_airfoil_unique.X_C-x_c));
            
            closest=new_upper_airfoil_unique(sortIndex(1),:);
            secound_closest=new_upper_airfoil_unique(sortIndex(2),:);
            
            average=interp1([closest.X_C,secound_closest.X_C], [closest.CP,secound_closest.CP], x_c, 'linear','extrap');
            
            %% saved with closest Ma and Cl as they are the same as in secound closest
            Temp=[Temp;{x_c, average, closest.Ma, closest.CL_res}];
        else
            Temp=[Temp;FIND];
        end  
        
    end
    %% Interpolation
    cp_inter=IDW(Temp.Ma, Temp.CL_res, Temp.CP, Ma, C_l, -2, 'ng', 4);
    
    Temp(:,:)=[];
    %% Result
    DV(i,:)=[x_c,cp_inter];
    i=i+1; 
end

9 个评论
显示 7更早的评论隐藏 7更早的评论

Gerrit 2020-1-9

在 MATLAB Online 中打开

new_upper_airfoil.mat

All standard function except IDW. Its an interpolation method.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%% INVERSE DISTANCE WEIGHT %%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function[Vint]=IDW(xc,yc,vc,x,y,e,r1,r2)
%%%%%%%%%%%%%%%%%
%%% INPUTS
%xc = stations x coordinates (columns) [vector]
%yc = stations y coordinates (rows) [vector]
%vc = variable values on the point [xc yc]
%x = interpolation points  x coordinates [vector]
%y = interpolation points y coordinates [vector]
%e = distance weight
%r1 --- 'fr' = fixed radius ;  'ng' = neighbours
%r2 --- radius lenght if r1 == 'fr' / number of neighbours if  r1 =='ng'
%%% OUTPUTS
%Vint --- Matrix [length(y),length(x)] with interpolated  variable values
%%% EXAMPLES
%%% --> V_spa=IDW(x1,y1,v1,x,y,-2,'ng',length(x1));
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   Simone Fatichi -- simonef@dicea.unifi.it
%   Copyright 2009
%   $Date: 2009/06/19 $ 
%   $Updated 2012/02/24 $ 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Vint=zeros(length(y),length(x));
xc=reshape(xc,1,length(xc));
yc=reshape(yc,1,length(yc));
vc=reshape(vc,1,length(vc));
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if  strcmp(r1,'fr')                                                         % compareStrings, fixed radius oder nearest neighbour
    if  (r2<=0)                                                             % Kontrolle, ob Radius positiv 
        disp('Error: Radius must be positive')
        return
    end
    for i=1:length(x)                                                       
        for j=1:length(y)
            D=[]; V=[]; wV =[]; vcc=[];
            D= sqrt((x(i)-xc).^2 +(y(j)-yc).^2);                            % Euklidische Norm des Abstands zwischen bekannten und gesuchten Punkten
            if min(D)==0                                                    % Kontrolle, ob Abstand bei einem bekannten Punkt 0 wird
                disp('Error: One or more stations have the coordinates of an interpolation point')
                return
            end
            vcc=vc(D<r2);                                                    % wählt vc-Werte der Wertepaare aus, bei denen D<r2 gilt 
            D=D(D<r2);                                                      % Euklidische Norm, die kleiner als der gesetzte Radius ist; alle dazugehörigen Z-Werte  
            V = vcc.*(D.^e);                                                % Zähler 
            wV = D.^e;                                                      % Nenner 
            if isempty(D)                                                   % Wenn Abstandsvektor D leer ist V=NaN   
                V=NaN;                                                      
            else
                V=sum(V)/sum(wV);                                           % Enformel, Summe über alle Einträge
            end
            Vint(j,i)=V;                                                    % Ergebnis
        end
    end
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
else                                                                        % r1 ist gleich 'ng'
    if (r2 > length(vc)) || (r2<1)                                          % Kontrolle, ob so viele Nachbarn wie gesetzt verfügbar sind
        disp('Error: Number of neighbours not congruent with data')
        return
    end
    for i=1:length(x)           
        for j=1:length(y)
            D=[]; V=[]; wV =[];vcc=[];                                      % erstellen der Vektoren
            D= sqrt((x(i)-xc).^2 +(y(j)-yc).^2);                            % Berechnen der euklidischen Norm
            if min(D)==0                                                    % Kontrolle, ob gesuchter Punkt in der Menge enthalten 
                disp('Error: One or more stations have the coordinates of an interpolation point') 
                return
            end
            [D,I]=sort(D);                                                  % sortiert D nach aufsteigender Reihenfolge, I ist Reihenfolge von vorher                                                
            vcc=vc(I);                                                      % vcc ist der Vektor mit den Z- Werten der Nachbarn in der Reihenfolge wie D
            V = vcc(1:r2).*(D(1:r2).^e);                                    % vcc(1:r2): Vektor mit z-Werten der nächsten Nachbarn; V = Zähler                                
            wV = D(1:r2).^e;                                                % Nenner der Endformel
            V=sum(V)/sum(wV);                                               % Endformel, Summe über alle Einträge
            Vint(j,i)=V;                                                    % Ergebnis
        end
    end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
return

Gerrit 2020-1-10

编辑：Gerrit 2020-1-10

在 MATLAB Online 中打开

I use another function now, so its save that I use the left and right neighours and not extrapolate.

t0=tic();
load('new_upper_airfoil');
load('new_lower_airfoil');
%% Eingabe der Parameter
Ma=0.7;
C_l=0.65;
% x_delta=0.008;
cl_delta=0.02;
Ma_delta=0.02;
%% Interpolation of upper_airfoil
% Delta=2/height(originalDV);
DV=zeros(202,2);
Temp=array2table(zeros(0,4),'VariableNames',{'X_C','CP','Ma','CL_res'});
Temp_update=array2table(zeros(0,4),'VariableNames',{'X_C','CP','Ma','CL_res'});
i=1;
%% load only the points within a Mach Number and Lift coefficient delta
new_upper_airfoil=new_upper_airfoil(new_upper_airfoil.Ma < (Ma+Ma_delta) & new_upper_airfoil.Ma > (Ma-Ma_delta) & new_upper_airfoil.CL_res < (C_l+cl_delta) & new_upper_airfoil.CL_res > (C_l-cl_delta),:);
new_lower_airfoil=new_lower_airfoil(new_lower_airfoil.Ma < (Ma+Ma_delta) & new_lower_airfoil.Ma > (Ma-Ma_delta) & new_lower_airfoil.CL_res < (C_l+cl_delta) & new_lower_airfoil.CL_res > (C_l-cl_delta),:);
for x_c=0:0.01:1
    %% Vector with all unique Ids (muss nicht nochmal für Unterseite erstellt werden)
    DataId_unique=unique(new_upper_airfoil.AeroDataId); %könnte gespeichert werden für alle Daten
    t_start=tic();
    for a=1:1:size(DataId_unique)
        %% temp table with all the information for one ID (AeroData ID is no longer needed)
        new_upper_airfoil_unique=new_upper_airfoil(new_upper_airfoil.AeroDataId==DataId_unique(a,:),2:5); % könnte später gespeichert werden
        
        %% try to find x/c
        FIND=new_upper_airfoil_unique(abs(new_upper_airfoil_unique.X_C-x_c)<0.0005,:); 
        
        if(isempty(FIND))
            %% find left and right neighbour
            leftNeighbourIndex=nearestpoint(x_c,new_upper_airfoil_unique.X_C,'previous');
            rightNeighbourIndex=nearestpoint(x_c,new_upper_airfoil_unique.X_C,'next');
            
            if isnan(leftNeighbourIndex)
                leftNeighbourIndex=rightNeighbourIndex+1;
            end
            if isnan(rightNeighbourIndex)
                rightNeighbourIndex=leftNeighbourIndex-1;
            end
            %% interpolation of the coordinate
            average=interp1([new_upper_airfoil_unique.X_C(leftNeighbourIndex,:),new_upper_airfoil_unique.X_C(rightNeighbourIndex,:)], [new_upper_airfoil_unique.CP(leftNeighbourIndex,:),new_upper_airfoil_unique.CP(rightNeighbourIndex,:)], x_c, 'linear','extrap');
            new_upper_airfoil_unique.X_C(1,:)=x_c;
            new_upper_airfoil_unique.CP(1,:)=average; 
            
            Temp=[Temp;new_upper_airfoil_unique(1,:)];
        else
            FIND.X_C=x_c;
            Temp=[Temp;FIND];
        end  
        
    end
    t_stuetz=toc(t_start);
    fprintf('Dauer für Stützstellen: %6.5f s\n',t_stuetz);
     
    cp_inter=IDW(Temp.Ma, Temp.CL_res, Temp.CP, Ma, C_l, -2, 'ng', 4);
    Temp(:,:)=[];
    DV(i,:)=[x_c,cp_inter];
    i=i+1;
    
end