Заполнение матрицы для для l oop в зависимости от того, является ли индекс четным или нечетным

Question

Это часть проблемы PDE, которую я пытаюсь решить численно. Я знаю, что мои расчеты для записей верны, однако я изо всех сил пытаюсь правильно собрать матрицу в Matlab. Первые строки верны, однако две последние строки не верны, так как мне не хватает элементов a_ {11,10}, a_ {10,10} и a_ {10,11}. Я скучаю по ним, потому что мой l oop никогда не обращается к ним. Например, для a_ {10,10} требуется i = 10, что является четным числом, однако мое M = 11 в этом случае и я иду только от 1 до 9. Если я сделаю мой для l oop go из 1:11, я выхожу за пределы ошибки, так как я использую индексы i + 1 и i + 2.

То же самое касается вектора b. Как мне заполнить все мои b_i, если я иду только к 9, но вектор имеет длину 11 элементов?

Есть предложения о том, как это исправить? Я попытался сделать отдельные циклы for для нечетного и четного случая, но у меня все еще остается та же проблема пропущенных записей из-за использования индексов i + 1 и i + 2. Или я должен вручную ввести (жестко) эти значения за пределами поля l oop?

Вот мой l oop:

---

h = 1/10; 
x = 0:h:1;
M = length(x);  

% Assemble stiffness matrix A and load vector b.

A = zeros(M,M);
b = zeros(M,1);

for i = 1:M-2 
    if rem(i,2) ~= 0                   % if i is odd
        A(i,i) = A(i,i) + 1/(3*h);

        b(i) = b(i) + f(x(i)) * (-2)*h/3;

    elseif rem(i,2) == 0               % if i is even
        A(i,i) = A(i,i) + 3/(3*h);
        A(i+2,i) = A(i+2,i) - 8/(3*h);
        A(i,i+2) = A(i,i+2) - 8/(3*h);

        b(i) = b(i) + f(x(i)) * 8*h/3;
    end

    A(i+1,i) = A(i+1,i) + 2/(3*h);
    A(i,i+1) = A(i,i+1) - 3/(3*h);
end

% Enforce the boundary conditions

A(1,1) = 1.e4;
A(end,end) = 1.e4;

---

Answer 1

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

Посмотрите на следующую модифицированную версию вашего кода:

h = 1/10; 
x = 0:h:1;
M = length(x);  

% Assemble stiffness matrix A and load vector b.

%Initialize A to have two more rows and columns.
A = zeros(M+2);
b = zeros(M,1);

for i = 1:M 
    if rem(i,2) ~= 0                   % if i is odd
        A(i,i) = A(i,i) + 1/(3*h);

        b(i) = b(i) + f(x(i)) * (-2)*h/3;

    elseif rem(i,2) == 0               % if i is even
        A(i,i) = A(i,i) + 3/(3*h);
        A(i+2,i) = A(i+2,i) - 8/(3*h);
        A(i,i+2) = A(i,i+2) - 8/(3*h);

        b(i) = b(i) + f(x(i)) * 8*h/3;
    end

    A(i+1,i) = A(i+1,i) + 2/(3*h);
    A(i,i+1) = A(i,i+1) - 3/(3*h);
end

%Remove two rows and columns form matrix A.
A = A(1:end-2, 1:end-2);

% Enforce the boundary conditions

A(1,1) = 1.e4;
A(end,end) = 1.e4;


function z = f(t)
% Stab f functionn - f(x) return x
z = t;
end

Я надеюсь, что код ведет себя так, как вы ожидали (не очень понятно, какое значение вы ожидаете прочитать за пределами A) ...