Постановка задачи.

Имеем систему линейно–независимых функций F(t)вL²[0 T]. На основе этих функций построить базис Ф(t). Для этого используем метод ортогонализации Грама–Шмидта [1]:

Здесь (f(t), f(t))— скалярное произведение вL²[0 T]. Из выражения следует:

Тогда задача сводится к нахождению коэффициентов матрицы ортогонализации c.

Решение.

Если решать задачу “в лоб” (**), то реализация функции в MATLAB будет такой:

function c = ortoC(t)

% Возвращает коффициенты ортогонализации [c]

% функций f[i](t) на промежутке [t].

% Кол-во функций определяется размером вектора

времени t

%

% @param t – время [0:dt:T]

% @return c – коэфф. ортогонализации

%

l = length(t);

c = zeros(l);

disp('Hello...');

% Скалярное произведение (f[i], f[j])

iFF = intFF(t);

% Поиск кэффиц. c[i,j]

for i = 1:l

% disp(i);

for j = 1:i

if (i == j)

c(i,j) = 1;

else

for nu = j:i-1

d_nu = 0;

for pi =1:nu

for pj =1:nu

d_nu = d_nu+c(nu,pi)*c(nu,pj)*iFF(pi,pj);

end;

end;

p_nu = 0;

for pl = 1:nu

p_nu = p_nu+iFF(i,pl)*c(nu,pl);

end;

c(i,j) = c(i,j)+p_nu/d_nu*c(nu,j);

end

c(i,j) = -c(i,j);

end;

end;

end;

disp('Bye');

% end function

В этой реализации алгоритма (**) у нас цикл имеет вложение пятого порядка, т.е. сложность алгоритма P определяется примерно так: P(l) = (a * l)⁵, гдеaнекоторый коэффициент. Поэтому при незначительном увеличении числа базисных функцийlвремя вычислений растёт пропорционально пятой степени.

Такая реализация конечно мало пригодна, тогда имеет смысл преобразовать (**) в матричную форму (***):

Тогда, с учётом (***), реализация (*) будет такой:

function [c, norm] = ortoC(t)

% Возвращает коффициенты ортогонализации [c] и нормировки [norm]

% функций f[i](t) на промежутке [t].

% Кол-во функций определяется размером вектора времени t

%

% @param t – время [0:dt:T]

% @return c – коэфф. ортогонализации

% norm – коэфф. нормировки

%

l = length(t);

c = zeros(l);

norm = zeros(1,l);

p = zeros(l);

disp('Hello...');

% Скалярное произведение (f[i], f[j])

iFF = intFF(t);

% Поиск кэффиц. c[i,j]

for n=2:l

%disp(n);

c(n, n) = 1;

% скалярное произведение

norm(n-1) = 1 / (c(n-1,:)*iFF*c(n-1,:)');

p(n,1:n-1) = (c(1:n-1,:)*iFF(:,n))' .* norm(1:n-1);

% элементы матрицы ортогонализации

c(n,1:n-1) = -p(n,1:n-1)*c(1:n-1,1:n-1);

end;

% коэф. нормировки

norm(l)=1 / (c(l,:)*iFF*c(l,:)');

disp('Bye');

% end function

Для этой реализации сложность алгоритма пропорциональна второй степени (если заменить операции перемножения матриц поэлементными циклами). Кроме того, операция умножения матриц реализована в ядре MATLAB, и выполняется несравненно быстрее цикла.

Алгоритм (*) используется для решения задач в области теории управления — “Метод моментов решения интегральных уравнений 1-го рода” [2]. При решении этой задачи были получены следующие временные показатели для числа базисных функций l=30:

	Цикл	Матричные операции
Интерпретатор	62,45 сек.	0,55 сек.
MEX–функция	22,79 сек.	0,68 сек.