Matlab: adapting myODEfun for the 'Vectorized' propery of odeset

Question

I have a stiff system of coupled ODEs that I am feeding MATLAB's ode15s solver. It works well, but now I'm trying to optimize the speed of integration. I am modeling 5 different variables on N different spatial sites, giving 5N coupled equations. For the moment, N=20 and integration time is about 25s, but I would like to go to larger values of N.

I used the profiler to see that the vast majority of the time is spent evaluating myODEfun. I did my best to optimize the code, but that doesn't change the fact that there is quite a bit going on in the function and that it is being evaluated ~50,000 times. I read that using the 'Vectorized' property for the ODEfunction can reduce the number of evaluations needed.

But I don't quite understand what exactly it is that I need to change about my ODEfun to make it conform to what Matlab wants a 'vectorized' ODEfun to look like.

From the documentation I see that you can change the example Van der Pol system from its normal form:

function dydt = vdp1000(t,y)
dydt = [y(2); 1000*(1-y(1)^2)*y(2)-y(1)];

to the vectorized form:

function dydt = vdp1000(t,y)
dydt = [y(2,:); 1000*(1-y(1,:).^2).*y(2,:)-y(1,:)];

I don't understand exactly what this new matrix of y is supposed to represent, and how the size of the second dimension is even defined. I could almost live with just adding ",:" and not thinking about it, but I am running into problems because I am already doing some vector operations in my code.

Here is a simplified example of my current functions, not yet vectorized. It models 2 variables, making 2*N equations. Please don't try to make sense of the ODEs that are generated here: they don't. I am talking about the operations that are happening.

function dydt = exampleODEfun(t,y,N)

dydt = zeros(2*N,1);
dTdt = zeros(N,1);
dXdt = zeros(N,1);

T = y(1:N);
X = y(N+1:2*N);

a = [T(2:N).^2 T(2:N) ones(N-1,1)];
b = [3 5 -2];

dTdt(1:N) = 0;
dXdt(1) = 0;
dXdt(2:N) = a*b';

dydt(1:N) = dTdt;
dydt(N+1:2*N) = dXdt;

end

Obviously in the real function a lot more is going on, both for T and X. As you can see, dXdt(1) is a boundary condition and needs its own calculations.

Blindly passing odeset 'Vectorized','on' and just adding ",:" to all the indexes does not work. For example, what size do I need to initialize dTdt and dXdt to now? What do I do to the ones(N-1,1)? And what do I need to do to make (a*b') still make sense?

I am using Matlab R2006a.

Answer 1

From help odeset:

Vectorized - Vectorized ODE function [ on | {off} ]

Set this property 'on' if the ODE function F is coded so that 
F(t,[y1 y2 ...]) returns [F(t,y1) F(t,y2) ...].

For the van der Pol example:

without vectorization:

function dydt = vdp1000(t,y)             %// 'y' is passed as [y1 
                                         %//                   y2]

    dydt = [y(2);                        %// 'dydt' is computed as [y1´
            1000*(1-y(1)^2)*y(2)-y(1)]   %//                        y2´]
                                         %// where the ´ indicates d/dt

with vectorization:

function dydt = vdp1000(t,y)            %// 'y' is passed as [y11 y21 y31 ...
                                        %//                   y12 y22 y32 ...] 

    dydt = [y(2,:);                             %// 'dydt' is computed as 
            1000*(1-y(1,:).^2).*y(2,:)-y(1,:)]; %//   [y11´ y21´ y31´ ...
                                                %//    y12´ y22´ y32´ ...]

where the y1, y2, y3, etc. refer to different vectors y at the same time t that ode15s will use to compute the next step.

For your example, you have to take into account that the y you get passed is no longer a vector, but a matrix in which every column represents a different vector you need to compute the derivative of:

function dydt = exampleODEfun(t,y,N)

    %// Adjust sizes to meet size of y
    dydt = zeros(2*N, size(y,2));
    dTdt = zeros(N, size(y,2));
    dXdt = zeros(N, size(y,2));

    %// Adjust indices to extract proper vales of ALL vectors
    T = y(1:N,:);
    X = y(N+1:2*N,:);

    %// This sort of section is usually where all the "thought" goes into:
    %// you can't use a*b' anymore, so I sum over the third dimension of the 
    % 3D array I built from your original vector
    b  = [3 5 -2];
    ab = sum(cat(3, b(1)*T(2:N,:).^2, b(2)*T(2:N,:), b(3)*ones(N-1, size(y,2))), 3);

    %// and finish it off
    dTdt(1:N,:) = 0;
    dXdt(1,:) = 0;
    dXdt(2:N,:) = ab;

    dydt(1:N,:) = dTdt;
    dydt(N+1:2*N,:) = dXdt;

end