matlab - 根据索引是偶数还是奇数，使用 for 循环填充矩阵

Question

这是我试图以数字方式解决的 PDE 问题的一部分。我知道我对条目的计算是正确的，但是我正在努力在 Matlab 中正确组装矩阵。第一行是正确的，但最后两行不正确，因为我缺少元素 a_{11,10}、a_{10,10} 和 a_{10,11}。我错过了这些，因为我的循环从不访问它们。例如 a_{10,10} 需要 i = 10，这是一个偶数，但是在这种情况下我的 M = 11，并且 i 只从 1 变为 9。如果我让我的 for 循环从 1:11 开始，我因为我使用的是 i+1 和 i+2，所以出现越界错误。

向量 b 也是如此。如果 i 只达到 9 但向量的长度为 11 个元素，我该如何填写我的所有 b_i？

关于如何解决这个问题的任何建议？我已经尝试为奇数和偶数情况制作单独的 for 循环，但由于使用了 i+1 和 i+2，我仍然遇到相同的条目丢失问题。或者我是否必须在 foo 循环之外手动输入（硬编码）这​​些不雅内容？

这是我的循环：

---

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;

---

标签： matlab

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）...