¿Cómo mejorar la velocidad de este código con bucles dobles en Matlab? ¿Es posible?

votos
0

Estimada comunidad de desbordamiento de pila. ¿Es posible mejorar la velocidad de este código en Matlab? ¿Puedo utilizar la vectorización? Tenga en cuenta que tengo que usar en cada bucle de la función fsolve vpasolve o.

Aquí está el código que llama a la función con un bucle doble:

   for i=1:1000
        for j=1:1000
            SilA(i,j)=SolW(i,j);  
        end
    end 

Esta es la función que es llamada por el código de seguridad. Puedo vectorizar la entrada de w, xi y hacer que el código se ejecute más rápido?

function [ffw] = SolW(w,xi)

format long e;
z=0;mm=0.46;sop=80;
epit=0.1;voP=80.;
rho=2.1;aposC=0.1;aposD=0.1;
parEh=0.2*10^6;parEv=0.2*10^6;parGv=0.074*10^6;
parpos=0.35;hp=0.2;Ep=30*10^6;
parposV=0.20;ll=0.15;dd=2*ll;

 C11=(parEh*(parEv-parEh*parpos^2)/((1+parpos)*(parEv-parEv*parpos-2*parEh*parpos^2)))*(1+2*1i*aposD);
C13=(parEh*parEv*parpos/((parEv-parEv*parpos-2*parEh*parpos^2)))*(1+2*1i*aposD);
C33=((1-parpos)*parEv^2/(parEv-parEv*parpos-2*parEh*parpos^2))*(1+2*1i*aposD);
C44=parGv*(1+2*1i*aposD);
DD=(Ep*hp^3)/(12*(1-parposV^2));
a1=C44;
a2=C33;
c1=(C13+C44)*1i*xi;
c2=(C13+C44)*1i*xi;
b1=rho*w^2-C11*xi^2;
b2=rho*w^2-C44*xi^2;

 syms xr
rsol=vpasolve((a1*xr+b1).*(a2*xr+b2)-c1*c2*xr==0,xr);
rsol=eval(rsol);
r=[-sqrt(rsol)];
r1=r(1,:);
r2=r(2,:);


 Fdf1=sop*sqrt(2*pi/(1i*epit*xi))*exp(1i*(xi*voP+w)^2/(2*epit*xi));
BC11=C44*(r1-1i*xi*c2*r1/(a2*r1^2+b2));
BC21=(C13*1i*xi)-((C33*c2*r1^2)/(a2*r1^2+b2))+(DD*xi^4-mm*w^2+1i*aposC*w)*c2*r1/(a2*r1^2+b2);
BC12=C44*(r2-1i*xi*c2*r2/(a2*r2^2+b2));
BC22=(C13*1i*xi)-((C33*c2*r2^2)/(a2*r2^2+b2))+(DD*xi^4-mm*w^2+1i*aposC*w)*c2*r2/(a2*r2^2+b2);

 syms As1 As2;
   try
[Ass1,Ass2]=vpasolve(BC11*As1+BC12*As2==0,BC21*As1+BC22*As2+Fdf1==0,As1,As2);
       A1=eval(Ass1);
       A2=eval(Ass2);
   catch
       A1=0.0;
       A2=0.0;
   end


Bn1=-(c2*r1/(a2*r1^2+b2))*A1;
Bn2=-(c2*r2/(a2*r2^2+b2))*A2;
ffw=Bn1*exp(r1*z)+Bn2*exp(r2*z);

end
Publicado el 19/09/2018 a las 17:11
fuente por usuario
En otros idiomas...                            


2 respuestas

votos
1

Todo en su función, pero las llamadas a vpasolve, y try....pueden ser vectorizar.

En primer lugar, todo esto no depende de wo xi, por lo que podría ser calculado por una sola vez :

format long e;
z=0;mm=0.46;sop=80;
epit=0.1;voP=80.;
rho=2.1;aposC=0.1;aposD=0.1;
parEh=0.2*10^6;parEv=0.2*10^6;parGv=0.074*10^6;
parpos=0.35;hp=0.2;Ep=30*10^6;
parposV=0.20;ll=0.15;dd=2*ll;

C11=(parEh*(parEv-parEh*parpos^2)/((1+parpos)*(parEv-parEv*parpos-2*parEh*parpos^2)))*(1+2*1i*aposD);
C13=(parEh*parEv*parpos/((parEv-parEv*parpos-2*parEh*parpos^2)))*(1+2*1i*aposD);
C33=((1-parpos)*parEv^2/(parEv-parEv*parpos-2*parEh*parpos^2))*(1+2*1i*aposD);
C44=parGv*(1+2*1i*aposD);
DD=(Ep*hp^3)/(12*(1-parposV^2));
a1=C44;
a2=C33;

Por lo que sé, evales lento, y estoy bastante seguro de que usted no lo necesita:

rsol=eval(rsol);

Aquí está un ejemplo de vectorización. Primero debe generar todas las combinaciones variables empleando meshgrid, y luego usar el .de Matlab notado utilizar operaciones prudentes elementos:

[I, J] = meshgrid(1:1000, 1:1000)
c1=(C13+C44)*1i.*xi;
c2=(C13+C44)*1i.*xi;
b1=rho.*w.^2 - C11.*xi.^2;
b2=rho.*w.^2-C44.*xi.^2;

Sin embargo, usted no será capaz de vectorizar vpasolve y try....litteraly, sin cambiar a otra cosa. vpasolvees probable que el cuello de botella de ustedes cómputo ( verificar esto utilizando perfiles de MATLAB ), optimizando así como se propone más arriba, probablemente no reducir el tiempo de cálculo mucho más.

A continuación, usted tiene varias soluciones:

  • usar parforsi usted tiene acceso a ella para paralelizar sus cálculos, que dependiendo de su arquitectura, que podría dar un aumento de velocidad 4x o menos.
  • puede ser posible para linealizar las ecuaciones y resolver todas en una sola operación. De todos modos, el uso de un programa de solución lineal será probablemente mucho más rápido que usar vpasolve. Esto probablemente le dará la aceleración más rápida (factor de adivinanzas 10 -100?)
  • porque tiene los datos continuos, se podría reducir el número de pasos, si se atreve precisión perder.

Espero que esto ayude

Respondida el 20/09/2018 a las 11:28
fuente por usuario

votos
0

En el programa de todo lo anterior se puede vectorizado como @beesleep dicho anteriormente. Por ejemplo:

[I, J] = meshgrid(1:1000, 1:1000)
c1=(C13+C44)*1i.*xi;
c2=(C13+C44)*1i.*xi;
b1=rho.*w.^2 - C11.*xi.^2;
b2=rho.*w.^2-C44.*xi.^2;   

La parte vpasolve, es decir,

 syms xr
 rsol=vpasolve((a1*xr+b1).*(a2*xr+b2)-c1*c2*xr==0,xr);
 rsol=eval(rsol);
 r=[-sqrt(rsol)];
 r1=r(1,:);
 r2=r(2,:);

También se pueden vectorizar utilizando fsolve como se muestra aquí:

fun=@(Xr) (a1.*Xr+b1).*(a2.*Xr+b2)-c1.*c2.*Xr
x01=-ones(2*Nxi);
x02=ones(2*Nxi);
options.Algorithm= 'trust-region-reflective'
options.JacobPattern=speye(4*Nxi^2);
options.PrecondBandWidth=0;
[rsol1]=fsolve(fun,x01,options);
[rsol2]=fsolve(fun,x02,options);

La otra parte, es decir,

 syms As1 As2;
 try
 [Ass1,Ass2]=vpasolve(BC11*As1+BC12*As2==0,BC21*As1+BC22*As2+Fdf1==0,As1,As2);
      A1=eval(Ass1);
      A2=eval(Ass2);
  catch
      A1=0.0;
      A2=0.0;
  end

desde contiene ecuaciones lineales y en cada fila tiene dos ecuaciones dependientes pueden ser resueltos como se muestra aquí:

 funAB1=@(As1) BC11.*As1+BC12.*((-Fdf2-BC21.*As1)./BC22);

 x0AB=ones(2*Nxi)+1i*ones(2*Nxi);
 options.Algorithm= 'trust-region-reflective';
 options.JacobPattern=speye(4*Nxi^2);
 options.PrecondBandWidth=0;
 [A1]=fsolve(funAB1,x0AB,options);

 A2=((-Fdf2-BC21.*A1)./BC22); 

Sin embargo, ambos también se pueden resolver analíticamente, es decir,

AAr=a1.*a2;
BBr=a1.*b2+b1.*a2-c1.*c2;
CCr=b1.*b2;

Xr1=(-BBr+sqrt(BBr^.2-4.*AAr.*CCr))./(2.*AAr);
Xr2=(-BBr-sqrt(BBr^.2-4.*AAr.*CCr))./(2.*AAr);
r1=-sqrt(Xr1(:,:));
r2=-sqrt(Xr2(:,:));
Respondida el 26/09/2018 a las 17:34
fuente por usuario

Cookies help us deliver our services. By using our services, you agree to our use of cookies. Learn more