Octave/MATLAB Заполнение элементов матрицы

Polyphem

Подскажите, каким образом можно убрать циклы из кода?
Допустим решается простейшая 3-x мерная дискретная задача на кубе с регулярной сеткой.
[math]$u - \tau \Delta u = f$[/math]
Эта задача очевидным образом сводится к СЛУ
[math]$K u = f$[/math]
Сейчас я делаю что-то подобное:

dim = (nx + 1) * ny * nz;
K = speye(dim, dim);
rhs = zeros(dim, 1);

for l = 1:nz
for m = 1:ny
for n = 2:nx
i = n + (m - 1) * (nx + 1) + (l - 1) * (nx + 1) * ny;

if (m == 1)
K(i, i - 1) = kix;
K(i, i + 1) = kix;
K(i, i + nx + 1) = kiy;
K(i, i) = 1 + kii + kiy;
elseif (m == ny)
K(i, i - 1) = kix;
K(i, i + 1) = kix;
K(i, i - nx - 1) = kiy;
K(i, i) = 1 + kii + kiy;
else
K(i, i - 1) = kix;
K(i, i + 1) = kix;
K(i, i) = 1 + kii;
K(i, i + nx + 1) = kiy;
K(i, i - nx - 1) = kiy;
end

if (l == 1)
K(i, i) = K(i, i) + kiz;
K(i, i + (nx + 1) * ny) = kiz;
elseif (l == nz)
K(i, i) = K(i, i) + kiz;
K(i, i - (nx + 1) * ny) = kiz;
else
K(i, i + (nx + 1) * ny) = kiz;
K(i, i - (nx + 1) * ny) = kiz;
end

rhs(i) = f(n, m, l);
end
end
end

Где K и rhs - матрица и правая часть СЛУ соответсвенно.
Условия внутри вложенных циклов - простые проверки
на то, лежит ли точка на правой/левой, нижней/верхней
границах куба.
Хочестся вместо трех вложенных циклов использовать
операции с векторами и разреженными матрицами.
Мне думалось, что это нечто стандартное, однако
нагуглить не получилось. Буду благодарен за помощь.

Niklz

некогда разбираться в коде, но такие многомерные сложносочиненные массивы в Матлабе обычно легко конструируются из одномерных векторов и стандартных операций при помощи функции bsxfun: http://www.mathworks.com/help/matlab/ref/bsxfun.html

stm7518204

Для матриц твоего вида, немеренно рулит следующий способ разбиения на блоки, где есть только нули снизу и блоки, которые содержат диагонали и значашие элементы.
A (1:10,1:10) = 1;
B = Zeros(10);
После чего полную матрицу ты собираешь как
K = [A,B;B,A] - 20 на 20
или большей размености.
K = [A,B,B;B,A,B; B,B,A] 30 на 30
Итд.
Известно, что в матлабе собрать циклами матрицу 200 на 84000 (: долго, а быстро было собрать раздельно 6 матриц по 200 на 14000 и потом воспользоваться одно операцией объединения. Выигрыш был очень большой. По крайней мере у меня в 2006 году.

Lene81

То, что я вижу у тебя — это построение матрицы в виде кронекеровского (поблочного) произведения матриц. В октаве/матлабе есть соотв. функция kron. Например, построение вклада от x переменной можно записать как kron(Tx, eye(n где Tx — матрица только в переменных x, а n — размер базиса
Далее, по одной переменной у тебя, как я вижу, ленточные матрицы. Их правильнее строить с помощью diag, которая может заполнять не только диагональ, но n-k над(под)диагонали. Например, если alpha — диагональ трехдиагональной матрицы, а beta_1 и beta_2 ее над- и поддиагонали, соответственно, то вся матрица собирается командой
A = diag(alpha) + diag(beta_1, 1) + diag(beta_2, -1)
Советую освоить эти команды (kron, diag, eye). Вначале можно писать "наивный" вариант с циклами, а потом проверять, как это сводится к матричным командам.

Polyphem

Благодарен всем отписавшимся. Буду разбираться.
Еще один маленький вопрос. Пусть функция f - функция аргументов (x, y, z, t).
Какой командой можно получить матрицу [math]$F(1:nx, 1:ny, 1:nz, 1:nt)$[/math] такую, что
[math]$F(i, j, k, l) = f(i \cdot hx, j \cdot hy, k \cdot hz, l \cdot ht)$[/math]
Сейчас я опять же использую цикл и команду feval.

makei

Если хватит памяти, то самый быстрый вариант такой:
[X, Y, Z, T] = ndgrid1:nx)*hx, (1:ny)*hy, (1:nz)*hz, (1:nt)*ht);
F = f(X, Y, Z, T);
где во второй строке ты оперируешь уже с матрицами, т.е. не забудь вместо * / ^ ставить соответственно *. /. ^.

makei

Лучше в начале посчитать всё без граничных условий, как я тебе написал до этого без циклов. А затем 4-мя строками кода поправить граничные условия.

Polyphem

Лучше в начале посчитать всё без граничных условий, как я тебе написал до этого без циклов
Не могу понять, что ты имеешь в виду? Не вижу твоего поста.
С матрицами еще предстоит разбираться, однако, воспользовавшись вашими советами/ссылками
в других местах программы удалось снизить время расчета с ~1.5 часов до ~20 минут. Продолжаю
дальнейшую оптимизацию :)
Оставить комментарий
Имя или ник:
Комментарий: