Как улучшить скорость этого кода с помощью двойных циклов в Matlab?Является ли это возможным?

Question

Уважаемое сообщество переполнения стека.Можно ли улучшить скорость этого кода в Matlab?Могу ли я использовать векторизацию?Обратите внимание, что я должен использовать в каждом цикле функцию «vpasolve» или «fsolve».

Вот код, который вызывает функцию с двойным циклом:

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

Вот функция, которая вызывается вышеуказанным кодом.Могу ли я векторизовать ввод w, xi и заставить код работать быстрее?

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

Answer 1

В вышеприведенной программе все может быть векторизовано, как сказано выше @beesleep.Например:

[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;   

Часть vpasolve, то есть

 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,:);

, также может быть векторизована с помощью fsolve, как показано здесь:

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);

другая часть, т. е.

 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

, поскольку содержит линейные уравнения и в каждой строке есть два зависимых уравнения, которые могут быть решены, как показано здесь:

 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); 

Однако оба также могут бытьрешено аналитически, т.е.

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(:,:));

Answer 2

Все в вашей функции, кроме вызовов vpasolve и try...., может быть векторизованным.

Во-первых, все это не зависит от 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;

Из того, что я знаю, eval медленный, и я уверен, что вам это не нужно:

rsol=eval(rsol);

Вот пример векторизации.Сначала вы должны сгенерировать все комбинации индексов, используя meshgrid, а затем использовать ., чтобы заметить matlab для использования поэлементных операций:

[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;

Но вы не сможете векторизовать vpasolve и try.... буквально, не меняя его на что-то другое.vpasolve вероятно является узким местом ваших вычислений ( проверьте это с помощью matlab profiler ), поэтому оптимизация, как предложено выше, вероятно, не сильно сократит ваше время вычислений.

Тогда у вас есть несколько решений:

• используйте parfor, если у вас есть доступ к нему для распараллеливания вычислений, что в зависимости от вашей архитектуры может дать вам 4-кратное ускорение или около того.
• может быть возможно линеаризовать ваши уравнения и решить их все за одну операцию.В любом случае, использование линейного решателя будет, вероятно, намного быстрее, чем использование vpasolve. Это, вероятно, даст вам самое быстрое ускорение (угадав коэффициент 10 -100?)
• , поскольку у вас есть непрерывные данные, вы можете уменьшить количество шагов, если осмелились потерять точность.

Надеюсь, это поможет