2015-09-21 2 views
0

Я пытаюсь использовать meshgrid в Matlab вместе с Chebfun, чтобы избавиться от двойных циклов. Я сначала определить квази-матрицу N функций,Векторизация с помощью Matlab Meshgrid в Chebfun

%Define functions of type Chebfun 
N = 10; %number of functions 
x = chebfun('x', [0 8]); %Domain 
psi = []; 
for i = 1:N 
psi = [psi sin(i.*pi.*x./8)]; 
end 

Расчет выборки будет вычислить двойную сумму $ \ sum_ {I, J = 1}^10 фунтов на квадратный дюйм (:, я). * Пси (:, J) $. Я могу добиться этого с помощью двух для петель в Matlab,

h = 0; 
for i = 1:N 
for j = 1:N 
h = h + psi(:,i).*psi(:,j); 
end 
end 

Затем я попытался использовать meshgrid векторизовать следующим образом:

[i j] = meshgrid(1:N,1:N); 
h = psi(:,i).*psi(:,j); 

я получаю ошибку «индекс столбца должен быть вектором целые». Как я могу решить эту проблему, чтобы избавиться от двойных циклов и сделать код более эффективным?

ответ

2

BTW, Chebfun не входит в состав родного MATLAB, и вам необходимо скачать его, чтобы запустить код: http://www.chebfun.org/. Однако это не должно влиять на то, как я отвечаю на ваш вопрос.

В принципе, psi - это колоночная матрица N, и это ваше желание складывать продукты из всех комбинаций пар столбцов в psi. У вас есть правильная идея с meshgrid, но вместо этого вы должны развернуть 2D-матрицу координат как для i, так и для j, чтобы они были единичными векторами. Затем вы должны использовать это и создать две матрицы столбцов N^2, которые соответствуют таким столбцам, где каждый столбец соответствует номерам точных столбцов, указанным в i и j, отобранных из psi. Затем вы выполняете умножение по элементам между этими двумя матрицами и суммируете все столбцы для каждой строки. BTW, я собираюсь использовать ii и jj в качестве переменных с выхода meshgrid вместо i и j. Эти переменные зарезервированы для комплексного числа в MATLAB, и я не хочу заслонять их непреднамеренно.

Что-то вроде этого:

%// Your code 
N = 10; %number of functions 
x = chebfun('x', [0 8]); %Domain 
psi = []; 
for i = 1:N 
psi = [psi sin(i.*pi.*x./8)]; 
end 

%// New code 
[ii,jj] = meshgrid(1:N, 1:N); 

%// Create two matrices and sum 
matrixA = psi(:, ii(:)); 
matrixB = psi(:, jj(:)); 
h = sum(matrixA.*matrixB, 2); 

Если вы хотите избавиться от временных переменных, вы можете это сделать в одном операторе после вызова meshgrid:

h = sum(psi(:, ii(:)).*psi(:, jj(:)), 2); 

Я не» t установить Chebfun, но мы можем проверить, что это вычисляет то, что нам нужно с помощью простого примера:

rng(123); 
N = 10; 
psi = randi(20, N, N); 

Выполнение этого кода с выше более эффективное решение дает нам:

>> h 

h = 

     8100 
     17161 
     10816 
     12100 
     14641 
     9216 
     10000 
     8649 
     9025 
     11664 

Кроме того, запуск выше двойной код for петли также дает нам:

>> h 

h = 

     8100 
     17161 
     10816 
     12100 
     14641 
     9216 
     10000 
     8649 
     9025 
     11664 

Если вы хотите быть абсолютно уверены, , мы можем иметь оба кода с выводами в виде отдельных переменных, а затем проверить, равны ли они:

%// Setup 
rng(123); 
N = 10; 
psi = randi(20, N, N); 

%// Old code 
h = 0; 
for i = 1:N 
    for j = 1:N 
     h = h + psi(:,i).*psi(:,j); 
    end 
end 

%// New code 
[ii,jj] = meshgrid(1:N, 1:N); 
hnew = sum(psi(:, ii(:)).*psi(:, jj(:)), 2); 

%// Check for equality 
eql = isequal(h, hnew); 

eql проверяет, равны ли обе переменные, и мы получаем их как таковые:

>> eql 

eql = 

    1 
Смежные вопросы