ODE15s running slowly for ODE with long expressions

조회 수: 2 (최근 30일)
Chen Wang
Chen Wang 2020년 1월 10일
댓글: darova 2020년 1월 10일
Dear Matlab experts,
I have a problem of ODE15s running slow. I am solving a differential equation in my research. This is the function of my equation:
function dp_dy=odefun(y,p)
global m m_ N f epsilon k alphaI AI t c
U=y+m_^2*epsilon^2/N*(cos(t.*(y-N))-1)./(y-N).^2;
U_y=1+m_^2*epsilon^2/N*( -t*(y-N).*sin(t*(y-N))-2*cos(t*(y-N))+2 )./(y-N).^3;
U_yy=m_^2*epsilon^2/N*( 4*t*(y-N).*sin(t*(y-N))-t^2*(y-N).^2.*cos(t*(y-N))+6*cos(t*(y-N))-6 )./(y-N).^4;
W= 2*real( conj(1i*m_^2*AI/2*( log(abs(y-N))-cosint(abs(y-N)*t)+1i*sinint((y-N)*t)-1i/2*pi+1 )+1i*AI*(N*alphaI-f)/(N^2-f*(f-1))).*(-1i*m_*N*AI/2*( (-1i*t*(y-N)-1).*exp(-1i*(y-N)*t)+1 )./(y-N).^2)/N^2 ...
+m_/N^3*abs(-1i*m_*N*AI/2*( exp(-1i*(y-N)*t)-1 )./(y-N)).^2 );
W_y= 2*real( conj(1i*m_^2*AI/2*( log(abs(y-N))-cosint(abs(y-N)*t)+1i*sinint((y-N)*t)-1i/2*pi+1 )+1i*AI*(N*alphaI-f)/(N^2-f*(f-1))).*(-1i*m_*N*AI/2*( exp(-1i*(y-N)*t).*(-t^2*(y-N).^2+2i*t*(y-N)+2)-2 )./(y-N).^3)/N^2+ ...
m_/N^3*(-1i*m_*N*AI/2*( exp(-1i*(y-N)*t)-1 )./(y-N)).*conj(-1i*m_*N*AI/2*( (-1i*t*(y-N)-1).*exp(-1i*(y-N)*t)+1 )./(y-N).^2)+ ...
2*m_/N^3*(-1i*m_*N*AI/2*( (-1i*t*(y-N)-1).*exp(-1i*(y-N)*t)+1 )./(y-N).^2).*conj(-1i*m_*N*AI/2*( exp(-1i*(y-N)*t)-1 )./(y-N)) );
G_y= 2*real( -1i*m_*conj(m_*AI/2*( exp(-1i*(y-N)*t)-1 )./(y-N)).*(m_*AI/2*( (-1i*t*(y-N)-1).*exp(-1i*(y-N)*t)+1 )./(y-N).^2)...
-conj(1i*m_^2*AI/2*( log(abs(y-N))-cosint(abs(y-N)*t)+1i*sinint((y-N)*t)-1i/2*pi+1 )+1i*AI*(N*alphaI-f)/(N^2-f*(f-1))).*(m_*AI/2*( exp(-1i*(y-N)*t).*(-t^2*(y-N).^2+2i*t*(y-N)+2)-2 )./(y-N).^3) );
h= -(2*(k*U-k*c+m*W).*(k*U_y+m*W_y)+f*U_yy)./((k*U-k*c+m*W).^2-f*(f-U_y))+m*( ( (k*U-k*c+m*W).^2-N^2-(k*U-k*c+m*W).^2).*W_y-1i*G_y*(k*U-k*c+m*W) )./((k*U-k*c+m*W).^2-N^2)./(k*U-k*c+m*W);
l2=k*f*(2*(k*U-k*c+m*W).*(k*U_y+m*W_y)+f*U_yy)./(k*U-k*c+m*W)./((k*U-k*c+m*W).^2-f*(f-U_y))-k^2 ...
+m./((k*U-k*c+m*W).^2-N^2).*( -m*((k*U-k*c+m*W).^2-f*(f-U_y))+k*f*((k*U-k*c+m*W).*W_y+1i*G_y)./(k*U-k*c+m*W) ) ;
dp_dy=zeros(2,1);
dp_dy(1)=p(2);
dp_dy(2)=-h.*p(2)-l2.*p(1);
end
and this is the command to solve the equation by calling the function
m_=0.5;m=0.5;N=4/3;epsilon=0.01;k=1;alphaI=-0.047;AI=-0.01+0.002;t=30;c=1.42+0.2i;
options = odeset('RelTol',1e-3,'AbsTol',1e-3);
[yout, pout]=ode15s( @odefun,[-10 3], [1;sqrt(k^2+m^2)]/1000,options );
It takes a long time (more than 100 seconds) on my computer to compute the results (I may need even higher accuracy and many more computations). I can understand that computing the ode function odefun takes some time, but these are all elemetary computations which I suppose should not be a problem for advanced software like Matlab. Moreover, whenever I paused the computation to see what was happening, it always paused in files like 'sym', 'char', 'size', 'digits' which are Matlab functions and had nothing to do with my own ode. Is it that Matlab was trying to do some symbolic compuation to my complicated odefun that takes so much time? If this is the case, then I only need numerical results and don't need symbolic solution. How to make it run faster?
P.S., I can say that the stiffness of the ode is not what makes it slow. In fact, I know the effect of W W_y and G_y in odefun are relatively weak, and if I replace the long expressions of W W_y and G_y with
W=0;W_y=0;G_y=0;
then the solution will be similar but the computation completes within 1 second. So it seems awkward to me that Matlab spent much more time in evaluating the function of the ode than solving the ode. I wonder if there is anyway to circumvent this difficulty.
Thank you very much for your time and your help!
Sincerely,
Chen

채택된 답변

darova
darova 2020년 1월 10일
편집: darova 2020년 1월 10일
  • Moreover, whenever I paused the computation to see what was happening, it always paused in files like 'sym', 'char', 'size', 'digits' which are Matlab functions and had nothing to do with my own ode
Indeed MATLAB uses some symbolic calculations in functions sinint and cosint. So to make your computations faster you can build your own functions. According to MATLAB documentation:
So sinint and cosint can be replaced with:
gamma = double( vpa('eulergamma') );
function res = ci(x)
t1 = linspace(eps,x,100);
res = gamma + log(x) + trapz(t1,(cos(t1)-1)./t1);
end
function res = si(x)
t1 = linspace(eps,x,100);
res = trapz(t1,sin(t1)./t1);
end
Can't see the difference (built-in functions and si/ci):
Edit: function names corrected
  댓글 수: 2
Chen Wang
Chen Wang 2020년 1월 10일
Thank you very much, Darova. The sinint and cosint were indeed what made the computation slow. I used your function instead and it ran much faster. I had a larger domain and higher demand of accuracy for ci and si functions, so I may need to increase the grid number. I think I can work it out. Thank you!
darova
darova 2020년 1월 10일
My plesure!

댓글을 달려면 로그인하십시오.

추가 답변 (0개)

카테고리

Help CenterFile Exchange에서 Symbolic Math Toolbox에 대해 자세히 알아보기

태그

Community Treasure Hunt

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

Start Hunting!

Translated by