Asked by Arsalan
on 17 Nov 2012

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

Answer by Jan Simon
on 18 Nov 2012

Edited by Jan Simon
on 25 Nov 2012

Accepted answer

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.

Show 8 older comments

Arsalan
on 26 Nov 2012

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

Opportunities for recent engineering grads.

## 2 Comments

## Jan Simon (view profile)

Direct link to this comment:http://www.mathworks.com/matlabcentral/answers/54009#comment_111752

What is your question? Is your problem stiff?

## Arsalan (view profile)

Direct link to this comment:http://www.mathworks.com/matlabcentral/answers/54009#comment_111815

Yes Jan, its a stiff ode.