Attempt to vectorize a Runge-Kutta algorith

3 次查看（过去 30 天）
Geoff
Geoff 2012-6-5
Good morning everyone !
I'm currently playing with a code which consists of a scores of Matlab functions, and one of which relies on a Runge-Kutta algorith. I'd like to improve the length of time required to implement it, namely by vectorizing the already-existing iterative one.
I have been carrying out several attempts for more than a week, and unfortunately, the final result appears to be a dismal failure. If you catch a glimpse at the code written just after my comment, you could understand that my major issue remains the vectorization of 'statp'. Indeed, 'statp' takes part in the outter-loop, but all its parameters are implemented in the inner-loop.
Therefore, could you help me solving this problem by vectorizing this code ?
Here is the sketch of code : (tab, tbh_beta_abs, beta_abs, p0, phi, alfa dtanlam and curv have already been defined in other functions. pchip enables to carry out some interpolations)
gamma=1.4;
am=tab(1,:);
bm=tab(2,:);
len=sqrt((am-am(1)).^2+(bm-bm(1)).^2);
tub=len/(len(end));
statp(1)=10^5;
beta_abs_tub=pchip(tbh_beat_abs,beta_abs,tub);
f(1)=0;
for k=1:(length(tub)-1)
dh=tub(k+1)-tub(k);
dz=am(k+1)_am(k);
dr=bm(k+1)-bm(k);
dlen=len(k+1)-len(k);
for j=1:4
Z=flor(j/2);
hh=tub(k)+.5*dh*Z;
zz=am(k)+.5*dz*Z;
rr=bm(k)+.5*dr*Z;
statpp=statp(k)+.5*f(j)*Z;
p0_hhm=pchip(tub,p0,hh);
squaremak=((statpp/p0_hhm)^(-(gamma-1)/gamma)-1)*2/(gamma-1);
newphi=pchip(tub,phi,hh);
cosalfa=cos(alfa);
sinalfa=sin(alfa);
coslam=cos(newphi-alfa);
sinlam=sin(newphi-alfa);
tanlam=tan(newphi-alfa);
dtanlam_int=pchip(tub,dtanlam,hh);
curvature=pchip(tub,curv,hh);
beta_abs_hh=pchip(tbh_beta_abs,beta_abs,hh);
squaremaktan=(sqrt(squaremak)*sin(beta_abs_hh))^2;
squaremaketa=(sqrt(squaremak)*cos(newphi-alfa)*cos(beta_abs_hh))^2;
f(j+1)=gamma*statpp*dlen/(1-squaremaketa/coslam^2)* ((1-squaremaketa)*cosalfa*squaremaktan/rr+...+curvature/(coslam^3)); %I have shortened the formula
end
statp((k+1)=statp(k)+1/6*(f(2)+2*(f(3)+f(4))+f(5));
end
Thank you for your help !
  2 个评论
Oleg Komarov
Oleg Komarov 2012-6-5
Please provide some example inputs so that we can run and profile it for you.
Geoff
Geoff 2012-6-5
I don't remember asking this question! ;-)

请先登录，再进行评论。

请先登录，再回答此问题。

回答（0 个）

类别

在 Help Center 和 File Exchange 中查找有关 Loops and Conditional Statements 的更多信息

标签

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by