Vectorisation d'une EDO dans Octave / Matlab

Question

Si j'ai une ode et que je l'ai écrite de deux façons, comme ici :

function re=rabdab()
x=linspace(0,2000,2000)';
tic;
[T,Y] = ode45(@fun,[x],[0 1 1]);
[T,Y] = ode45(@fun,[x],[0 1 1]);
[T,Y] = ode45(@fun,[x],[0 1 1]);
toc;

tic;
[A,B] = ode45(@fun2,[x],[0 1 1]);
[A,B] = ode45(@fun2,[x],[0 1 1]);
[A,B] = ode45(@fun2,[x],[0 1 1]);
toc;

function dy = fun(t,y)
dy = zeros(3,1);    % a column vector
dy = [y(2) * y(3);...
-y(1) * y(3);...
-0.51 * y(1) * y(2);];

function dy = fun2(t,y)
dy = zeros(3,1);    % a column vector
dy(1) = y(2) * y(3);
dy(2) = -y(1) * y(3);
dy(3) = -0.51 * y(1) * y(2);

Il n'y a pratiquement pas de différence de temps. L'un prend autant de temps que l'autre. Mais je pensais que fun est la version vectorisée de fun2 . Ou est-ce que je me trompe ? Le but est d'accélérer un peu mon code. L'exemple est tiré de la page web de matlab. Je crois que je n'ai pas bien compris ce que signifie "vectorisé". Si c'est déjà vectorisé, à quoi ressemblerait un code non vectorisé ?

Answer 1

Vectorisation est un concept étroitement lié à celui de la programmation fonctionnelle . En MATLAB, il s'agit d'effectuer des opérations sur des tableaux (vecteurs ou matrices) sans écrire implicitement une boucle.

Par exemple, si vous deviez calculer la fonction f(x) = 2x pour tout entier x entre 1 et 100, vous pouvez écrire :

for x = 1:100
    f(x) = 2 * x;
end

qui n'est pas un code vectorisé. La version vectorisée l'est :

x = 1:100; %// Declare a vector of integer values from 1 to 100
f = 2 * x; %// Vectorized operation "*"

ou même plus court :

f = 2 * (1:100);

MATLAB étant un langage interprété, il est évident que l'interpréteur traduit cela en une sorte de boucle "sous le capot", mais il est optimisé et généralement beaucoup plus rapide que l'interprétation réelle d'une boucle (lisez cette question pour référence). Enfin, en quelque sorte, il en a été ainsi jusqu'aux versions récentes de MATLAB, dans lesquelles Accélération JIT a été intégré (lire ici ).

Revenons maintenant à votre code : vous avez ici deux versions vectorisées de votre code : l'une qui concatène verticalement trois valeurs et l'autre qui assigne directement ces valeurs dans un vecteur de colonnes. C'est fondamentalement la même chose. Si vous l'aviez fait avec une boucle for explicite, cela n'aurait pas été "vectorisé". En ce qui concerne le gain de performance réel de la "vectorisation" d'une boucle (c'est-à-dire la conversion d'une boucle for en une opération vectorisée), cela dépend de la rapidité de la boucle for en raison de l'accélération JIT en premier lieu.

Il ne semble pas qu'il y ait grand-chose à faire pour accélérer votre code. Vos fonctions sont assez basiques, donc cela se résume à l'implémentation interne de ode45 que vous ne pouvez pas modifier.

Si vous souhaitez en savoir plus sur la vectorisation et l'écriture de code MATLAB plus rapide en général, voici un article intéressant : Loren sur l'art de MATLAB : " Accélérer les applications MATLAB " .

Bon codage !

Answer 2

La réponse précédente est correcte d'un point de vue général, vectorisé dans le contexte des ODE signifie quelque chose de plus spécifique. En bref, une fonction f(t,y) est vectorisé si f(t,[y1 y2 ...]) retours [f(t,y1) f(t,y2) ...] , pour y1,y2 vecteurs de colonnes. Selon la documentation [1], "cela permet au solveur de réduire le nombre d'évaluations de fonctions nécessaires pour calculer toutes les colonnes de la matrice jacobienne".

Fonctions fun3 et fun4 ci-dessous sont correctement vectorisés au sens de l'EDO :

function re=rabdab()
x=linspace(0,20000,20000)';
opts=odeset('Vectorized','on');

tic;
[T,Y] = ode45(@fun,[x],[0 1 1]);
[T,Y] = ode45(@fun,[x],[0 1 1]);
[T,Y] = ode45(@fun,[x],[0 1 1]);
toc;

tic;
[A,B] = ode45(@fun2,[x],[0 1 1]);
[A,B] = ode45(@fun2,[x],[0 1 1]);
[A,B] = ode45(@fun2,[x],[0 1 1]);
toc;

tic;
[A,B] = ode45(@fun3,[x],[0 1 1],opts);
[A,B] = ode45(@fun3,[x],[0 1 1],opts);
[A,B] = ode45(@fun3,[x],[0 1 1],opts);
toc;

tic;
[A,B] = ode45(@fun4,[x],[0 1 1],opts);
[A,B] = ode45(@fun4,[x],[0 1 1],opts);
[A,B] = ode45(@fun4,[x],[0 1 1],opts);
toc;

function dy = fun(t,y)
dy = zeros(3,1);    % a column vector
dy = [y(2) * y(3);...
-y(1) * y(3);...
-0.51 * y(1) * y(2);];

function dy = fun2(t,y)
dy = zeros(3,1);    % a column vector
dy(1) = y(2) * y(3);
dy(2) = -y(1) * y(3);
dy(3) = -0.51 * y(1) * y(2);

function dy = fun3(t,y)
dy = zeros(size(y)); % a matrix with arbitrarily many columns, rather than a column vector
dy = [y(2,:) .* y(3,:);...
-y(1,:) .* y(3,:);...
-0.51 .* y(1,:) .* y(2,:);];

function dy = fun4(t,y)
dy = [y(2,:) .* y(3,:);... % same as fun3()
-y(1,:) .* y(3,:);...
-0.51 .* y(1,:) .* y(2,:);];

(Remarque annexe : si l'on omet l'allocation de mémoire inutile à l'aide de zeros laisse fun4 fonctionnent légèrement plus vite que les fun3 .)

• Q : Qu'en est-il de la vectorisation par rapport au premier argument ?
• R : Pour les solveurs ODE, la fonction ODE est vectorisée uniquement par rapport au deuxième argument. Le solveur de problème de valeur limite bvp4c nécessite cependant une vectorisation en ce qui concerne les premier et deuxième arguments. [1]

La documentation officielle [1] fournit des détails supplémentaires sur la vectorisation spécifique aux ODE (voir la section "Description des propriétés du jacobien").

[1] https://www.mathworks.com/help/releases/R2015b/matlab/ref/odeset.html