Speeding Up the ODE solver algorithm for a given time dependent input

Hi,

I am currently looking for a way to modify my code such that it will require shorter computation time, for solving a system of differential equations. The part of my code which is causing the most computation time in brief is the following ODE solver:

the ODE solver is called in my main program as

options=odeset('AbsTol',10e-6,'RelTol',10e-5,'Stats','on'); 
[T_, Y_]=ode45(@yourFcn, tspan,Y0, options, freq, b, amp);

where Y0 is a vector with initial conditions, freq, b, amp are constants passed to my function. tspan is a vector of X elements where in my case I like it to be approximately 50000 elements. and I have evaluated that I need the corresponding Absolute and relative error tolerances for my system.

function dy = yourFcn(t, y, b, freq, amp);
J=(amp/2)*(sin(2*pi*freq*t);
J=J+b;

dy = zeros(3,1); 
dy(1) = (A*(y(3))/(B*y(1)))*y(1)+ C*y(3)^2;  
dy(2) = (B*(y(3))/(A*y(1))-1);               
dy(3) = J/C*y(3) - D*(y(3))/(A*y(1))*y(1);

where A,B,C,D are all constants, and J is my time dependent variable in my system of differential equations which will also have a length of approximately 50000 elements as same as tspan.

Looking at the profile function, as expected most of time has been taking by calling the function yourFcn and the algorithem used for ODE45.

I have so far tried to evaluate the vector J in my main program, evaluating each element of J over a short tspan of a:b, and capturing the steady state solution of my differential equation as my final result for that particular element of J and moving to next point......

I have also tried the following method by defining J in my main program again and calling my ODE solver as follow

[T_, Y_]=ode45(@(tspan, y) yourFcn(tspan, y, tspan, J) ,tspan, Y_0, options);

where in my Function i made the following changes

function dy = yourFcn(t, y, t1, J1);
dy = zeros(3,1); 
J=interp1(t1, J1, t);
.
.
.

but in all case the first method showed the best performance but still extremely slow for me.

My computer specs are: Intel Quad Core i7 @ 3.7Ghz 4Gb Ram

Many Thanks, Arsalan

[Jan]

What is your question? Is your problem stiff?

[Arsalan]

Yes Jan, its a stiff ode.

Answer 1

[Jan]

ODE45 is not sufficient for stiff ODEs. ODE15S or ODE23S is better.

[EDITED] If you have a mutli-core machine, you cannot simply distribute the integration by splitting the interval, because you need the final value of the previous value top start. But Matlab's integrators profit from vectorised integrands, such that t is a vector and y a matrix.

I cannot understand this sentence: "I have predefined my time tspan as a vector where dt or (delta_t) would be 1/(sampling_rate) where my sampling rate is around 150e9 and I have about 50000 elements in my tspan vector". If dt is 1/150e9 and you have 50000 elements, does this mean, that the simulation time is 6.66e-7 seconds? Although this might be the physical reality, it can cause unexpected rounding errors. I suggest to scale the time such that you get something near to 1.

Again I insist, that ODE45 is a bad choice for a stiff system. If the stiff solvers are slower, you should exhaustively check, if there is another bug. Using a non-stiff solver for a stiff problem is very likely a numerical mistake and I would not accept a paper or PhD, when this is done without really very good reasons and a detailed analysis of stability. computing time and the claim that the accuracy is "sufficient" are no hard scientific arguments.

Unfortunately Matlab's integrators handle tspan by evaluating the time to the specified elements of tspan. Smarter integrators keep their step size provided by the step size control and get the values to tspan by an interpolation using the already evaluated intermediate positions. Depending on how dynamic your system is, you can simulate this externally: Use 5000 instead of 50000 time steps in tspan and interpolate the output.

[EDITED 2]

A "vectorized odeFunc" reduces the computing time. From "doc ode23":

 Set the Vectorized property 'on' if the ODE function is coded so that odefun(T,[Y1,Y2 ...]) returns [odefun(T,Y1),odefun(T,Y2) ...].

Another improvement: Use the JPattern property in odeset() to exploit the sparsity pattern of the mass matrix (if there is one).

Set the MStateDependence flag to consider a weak or no dependency between y and the Mass-Matrix.

Providing the Jacobian explicitly can be much faster than letting the integrator evaluate it dynamically by finite differences.

qinterp1 can be improved, if you mean FEX: qinterp1.m. You could try http://www.mathworks.com/matlabcentral/fileexchange/25463-scaletime, but this is actually designed for larger problems. It seems promising to remove all unnecessary code from qinterp1 instead.

But finally there could be a design problem also: J=interp1(t1, J1, t) seems to be a linear interpolation inside the function to be integrated. As far as I can see, the linear interpolation is not differentiable. But Matlab's integrators required a smooth function. Otherwise the discontinuities in the derivatives squeeze the step size control out of the specifications. While you will still get a result, it has not been calculated in a scientifically clean way. Driving the integrator on a non-smooth function will reduce the accuracy of the result at first by producing bad derivatives and by accumulating more rounding errors due to choosing too small step sizes.

Nevertheless, integrating discontinuous functions is very common and I've seen this in published papers as well as in PhD theses. Although others do this frequently, it is a bad idea.

[Arsalan]

for this particular system that I am using, it seems as I am getting enough accuracy via ode45, but substituting the solver to either 15s or 23s, will lead to much higher computation times.

[Jan]

Then please explain again your question. Actually using ODE45 on a stiff ODE leads to tiny stepsizes, such that the result contains a large accumulated roundoff error.

[Arsalan]

I have predefined my time tspan as a vector where dt or (delta_t) would be 1/(sampling_rate) where my sampling rate is around 150e9 and I have about 50000 elements in my tspan vector. as you can see this by itself is increasing the computation time as many points is to be evaluated.

Speed is really important for me, as I will have to test and run the simulation over many times. thus I am looking for a way that I can speed up my work. It looks as Ode45 is giving me enough accuracy and speed compared to ode15s and 23s, even though its taking smaller step sizes.

Giving the methods that I have defined my problem do you think is there any way that I can rearrange the code such it will run faster. Or the only solution would be to use a more powerful computer or taking advantage of Parallel computing or Mex's.

Thanks, Arsalan

[Arsalan]

you are right it will equate to some fraction of second in order of 1e-7.

I'll try to compare the results I am acheiving from my ODE solver to an anlytical solution, to see if much errors are produced, I'll let you know how it goes on. Thanks

[Arsalan]

Hi, You were dead right about the use of ODE13s and ODE23s, they lead to more accurate results compared to analytical results. Now that I am using ODE23s, I am still having a problem with the computation speed as my varying input to the ODE solver "J" is a very long vector. What suggestions would you have for me to be able to speed up the computation Arsalan

[Jan]

Besides running a faster machine? Use the profiler to find the bottlenecks. Drink a cup of coffee (I'm not kidding - science involves patience). Is the function to be integrated vectorized?

[Arsalan]

Belive me I have been drinking allot of coffe for past few days. I've been using the profiler, at the start there was allot of time spent on calling the interp1 function which was positioned in the file containing the differential equations, I'm now using qinterp1 function, which has significantly helped the speed. the other bottelnecks are being caused by ode23s itself at functions

"feval" on lines (430, 442, 412, 443) and "odenumjac" at "feval" command (line 132)

But I am not very good at Matlab thus I dont really know how to change these algorithems to get higher performance.

I am currently using a fairly better machine with a dual core Xeon Processor. but the problem is that it seems as only one thread is being used at all times. I am currently trying to see if I can use the paralell toolbox to make full use of all threads available.

I didn't understand what you meant by "function to be integrated vectorized"?

Thanks

[Jan]

"A lot of coffee" does not solve the problem. The trick is to drink more coffee. ;-)

See [EDITED 2] in my answer.

[Arsalan]

hahaha, thanks for the tip, I'll start drinking more from now on. Thanks for the comments I'll try working on them, I'll let you know about the results

[Arsalan]

Hi Jan, Found the first problem which I cant fix, so I tired to vectorize "myFcn" as follow

making sure that initial conditions are = [x, y, z]

function dy = yourFcn(t, y, b, freq, amp);
J=(amp/2)*(sin(2*pi*freq*t);
J=J+b;
dy = zeros(3,1); 
  dy(1, :) = (A*(y(3, :))/(B*y(1, :)))*y(1, :)+ C*y(3, :)^2;  
  dy(2, :) = (B*(y(3, :))/(A*y(1, :))-1);               
  dy(3, :) = J/C*y(3, :) - D*(y(3, :))/(A*y(1, :))*y(1, :);

where the following error showed

Subscripted assignment dimension mismatch.

I tried the following way

function dy = yourFcn(t, y, b, freq, amp);
J=(amp/2)*(sin(2*pi*freq*t);
J=J+b;
dy = zeros(3, 1); 
  dy(:,1) = (A*(y(:,3))/(B*y(:,1)))*y(:,1)+ C*y(:,3)^2;  
  dy(:,2) = (B*(y(:,3))/(A*y(:,1))-1);               
  dy(:,3) = J/C*y(:,3) - D*(y(:,1))/(A*y(:,1))*y(:,1);

and the following error showed up

Attempted to access y(:,3); index out of bounds because size(y)=[3,1].

What do you think is wrong here. Thanks

[Jan]

You need the elementwise operations .*, ./ and .^, if you do not mean the matrix operations.

y(3, :) ^ 2

becomes:

y(3, :) .^ 2