How to speed up the 'dde23' solver?

Dear all,

I try to use the solver 'dde23' to simulate the time-dealy differential equation with single constant delay. But the running time of the solver is too long. Is there any way to speed up the solver without reducing the calculation accuracy?

The main part of my code is attached as follows:

clear all
clc

kc=0.884;  
lags=[1.2265834]; 
tspan=[0,2000];

options=ddeset('RelTol',1e-8,'AbsTol',1e-8,'MaxStep',0.05);
sol=dde23(@(t,y,Z)ddefun(t,y,Z,kc),lags,@(t)history(t),tspan,options);

[Jan]

Please show us the ddefun to give us a chance to see, if it can be improved.

Most of all omit the darn clear all , which is a waste of time: It removes all loaded functions from the memory and reloading them from disk wastes a lot of time. clc is not useful also, is it?

Of course there is no magic "FasterCalculations" flag settable by odeset.

[Shan Yin]

Dear Jan,

Thanks for your comment. The rest subfunctions are attached as follows. These functions are applied to simulate the bifurcating response near a degenerate Hopf bifurcation point. Can they be improved?

function h=history(t)
   w0=0.809535221224209;  
   r=0.00001;   
   q0=[1;0.809535221224209*1i;0.914698562098013+0.190646906623320e-1*1i];
   x1=2*real(q0(1)*exp(1i*w0*t))*r;
   x2=2*real(q0(2)*exp(1i*w0*t))*r;
   x3=2*real(q0(3)*exp(1i*w0*t))*r;
   h=[x1;x2;x3];
end


function dydt=ddefun(t,y,Z,kc) 
  n=3;
  e=0.3;
  b=0.4;
  J=1.2;
  g=9.8;
  l=0.5;
  r=0.05;
  gama=0.7;
  beta=0.49;
  a=r/l;
  x0=acos(beta/gama);
  z0=sqrt(((e*gama+g*gama/l-e*beta)*sqrt(gama^2-beta^2))/(beta*n^2*sqrt(gama^2-beta^2)+beta*a*n^2*gama));
  a0=(e+n^2*z0^2)/2;
  a1=n^2*z0;
  a2=n^2/2;
  a3=a*n^2*z0^2;
  a4=2*a*n^2*z0;
  a5=a*n^2;
  a6=e+g/l;
  a7=gama/J;
  a8=beta/J;
  
  dydt=[y(2);
       (a0+a1*y(3)+a2*(y(3))^2)*sin(2*y(1)+2*x0)+(a5*(y(3))^2+a4*y(3)+a3)*cos(y(1)+x0)-a6*sin(y(1)+x0)-b*y(2);
        a7*cos(y(1)+x0)-a8+kc*(y(1)-Z(1,1))];
end

[Jan]

|exp| is a very expensive operation (but this is not the reason for its name). So do not call it 3 times for the same input. You can replace:

   x1=2*real(q0(1)*exp(1i*w0*t))*r;
   x2=2*real(q0(2)*exp(1i*w0*t))*r;
   x3=2*real(q0(3)*exp(1i*w0*t))*r;
   h=[x1;x2;x3];

by

   h = 2 * real(q0 * exp(1i * w0 * t)) * r;

This is 1 call of exp() only and nicer.

As far as I can see, x0, z0, a0...a8 are constants. Especially

  z0=sqrt(((e*gama+g*gama/l-e*beta)*sqrt(gama^2-beta^2)) / ...
      (beta*n^2*sqrt(gama^2-beta^2)+beta*a*n^2*gama));

look expensive. Then define the variables once only and store them persistently:

  persistent z0
  if isempty(z0)
      z0 = sqrt(((e*gama+g*gama/l-e*beta)*sqrt(gama^2-beta^2)) / ...
                 (beta*n^2*sqrt(gama^2-beta^2)+beta*a*n^2*gama));
  end

[Steven Lord]

Then define the variables once only and store them persistently:

Or compute some of those parameter values outside the DDE function and pass them in as additional parameters.

[Jan]

@Shan Yin: Steven's suggestion is smarter than storing the parameters persistently inside the function. It is the nature of parameters to be defined separately. But bother methods should be a benefit for the run time also.

[Shan Yin]

@Jan @ Steven Lord. Thanks for your nice suggestions!