# ODE15s running slowly for ODE with long expressions

1 view (last 30 days)
Chen Wang on 10 Jan 2020
Commented: darova on 10 Jan 2020
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 on 10 Jan 2020
Edited: darova on 10 Jan 2020
• 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

Chen Wang on 10 Jan 2020
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 on 10 Jan 2020
My plesure!