Efficient ODE with function handles

Hi,

I am solving an ODE very often, so efficiency is important. The ODE I am solving calls other functions, because the ode is very general, but the details depend on the problem.

What I am doing now is this ( not working code, but the relevant snippets)

% just some parameters that are passed. This does not impact performance
params   = {1,2,3}; 
% function handle of function that changes the ode in question
h1  	 = @f_function;
% function handle of the ode 
%h_dydt   = @(t,y)f_dydt(t,y,params);   % aproach 1)
h_dydt   = @(t,y)f_dydt(t,y,params,h1); % aproach 2)
solution = ode45(h_dydt,...)
   

Where the ode function is

%function [dydt] = f_dydt(t,y,params,h)% aproach 1)
function [dydt] = f_dydt(t,y,params,h) % aproach 2)
% a = f_function(y,params); % aproach 1)
a =  h(y,params);             % aproach 2)
%the rest is not important,but could be 
dydt = y*t*a;
end

The problem is, that aproach 2) takes 60% more time than aproach 1), and I wonder if there is a way to have the same flexibility without the huge added cost. As it is now, I have to do aproach 1) and write a separate dydt function for each problem.

Thanks for looking into this!

Answer 1

[Bjorn Gustavsson]

If the function-handle-argument is too costly and you still want "some general flexibility" you might try to simply extend the ODE-function with a switch-selecting between the different drivers? Something like this:

%function [dydt] = f_dydt(t,X,params,h)% aproach 1)
function [dydt] = f_dydt(t,X,params,n_RHS) % aproach 1.5)
% a = f_function(x,y,params); % aproach 1)
switch n_RHS
  case 1
    a = h1(X,params);             % aproach 1.5)
  case 2
    a = h2(X,params);             % aproach 1.5)
  otherwise
    a = 12;
end
%the rest is not important,but could be 
dydt = y*t*a;
end

This should fall somewhere inbetween your two aproaches, and might be good enough as a solution to both issues. It will require some effort on writing the documentation of your function to provide "easy use".

HTH

[Stefan Weichert]

Hi,

thanks a lot for your suggestion. This is indeed a possibility.

Answer 2

[Steven Lord]

Does f_function expect to be called with two input arguments or three?

x = [1; 2; 3];
y = [4; 5; 6];
value1 = min([x, y], 4) % This works

value1 = 3×2

1 4 2 4 3 4

value2 = min(x, y, 4) % This does not work

Error using min
Dimension argument is not supported when two input arrays are provided.

How are the X input argument with which ode45 calls your function related to the two variables x and y whose definitions you've elided from your code?

I suspect you're comparing apples and oranges here, though with the limited information you've provided I can't be certain.

[Stefan Weichert]

I see I mixed x and y after reducing the complexity of the functions to declutter the code.

Fixed it now. Thanks.

But the problem with efficiency remains.