MATLAB求解二阶微分方程遇‘索引超出数组元素数目’错误求助

我用有限差分法在MATLAB中求解二阶微分方程，反复碰到Index exceeds the number of array elements错误。已经尝试调整findiff.m里Vb向量的长度，但问题依旧。需要定位错误根源并解决。

相关代码 #

trisys.m #

function X = trisys(A, D, C, B)
N = length(B);

% Forward elimination
for k = 2:N
    mult = A(k-1) / D(k-1);
    D(k) = D(k) - mult * C(k-1);
    B(k) = B(k) - mult * B(k-1);
end

% Back substitution
X(N) = B(N) / D(N);
for k = N-1:-1:1
    X(k) = (B(k) - C(k) * X(k+1)) / D(k);
end
end

findiff.m #

function F = findiff(p, q, r, a, b, alpha, beta, N)
T = zeros(1, N + 1);
X = zeros(1, N - 1);

Va = zeros(1, N - 2);
Vb = zeros(1, N - 1);
Vc = zeros(1, N - 2);
Vd = zeros(1, N - 1);
h = (b - a) / N;

Vt = linspace(a + h, b - h, N - 1);

Vb = -h^2 * feval(r, Vt);
Vb(1) = Vb(1) + (1 + h / 2 * feval(p, Vt(1))) * alpha;
Vb(N - 1) = Vb(N - 1) + (1 - h / 2 * feval(p, Vt(N - 1))) * beta;

Vd = 2 + h^2 * feval(q, Vt);

Vta = Vt(2:N - 1);
Va = -1 - h / 2 * feval(p, Vta);

X = trisys(Va, Vd, Vc, Vb);

T = [a, Vt, b];
X = [alpha, X, beta];
F = [T' X'];
end

p.m #

function result = p(t)
result = 2 * t / (1 + t^2);
end

q.m #

function result = q(t)
result = -2 / (1 + t^2);
end

r.m #

function result = r(t)
result = 1;
end

main.m #

a = 0;
b = 4;
alpha = 1.25;
beta = -0.95;
N = 100;

result = findiff(@p, @q, @r, a, b, alpha, beta, N);

disp(result);

错误根源定位 #

1. 元素运算缺失导致向量转矩阵：p.m和q.m中未使用MATLAB的元素-wise运算（.），当输入为向量时，t^2会执行矩阵乘法而非元素平方，返回值变为矩阵，赋值给Va、Vd后破坏了向量结构，后续trisys函数索引操作必然越界。
2. 上对角线向量未赋值：findiff.m中Vc（三对角矩阵的上对角线）始终为全零，既不符合有限差分的离散格式，也会导致方程求解逻辑错误。
3. 变量名混淆：trisys.m中用N表示未知量个数，与主函数的网格数N重名，容易引发逻辑混乱。

修复方案 #

1. 修正p.m和q.m的元素运算 #

将函数中的乘法、平方改为元素操作，确保向量输入时正确计算每个元素：

% p.m
function result = p(t)
result = 2 .* t ./ (1 + t.^2);
end

% q.m
function result = q(t)
result = -2 ./ (1 + t.^2);
end

2. 正确赋值findiff.m中的Vc向量 #

在Va赋值后添加上对角线元素的计算代码：

% findiff.m中，Va赋值下方添加
Vtc = Vt(1:N-2);
Vc = -1 + h / 2 * feval(p, Vtc);

3. （可选）修改trisys.m变量名避免混淆 #

将trisys中的N改为M，明确表示未知量个数：

function X = trisys(A, D, C, B)
M = length(B);

% Forward elimination
for k = 2:M
    mult = A(k-1) / D(k-1);
    D(k) = D(k) - mult * C(k-1);
    B(k) = B(k) - mult * B(k-1);
end

% Back substitution
X(M) = B(M) / D(M);
for k = M-1:-1:1
    X(k) = (B(k) - C(k) * X(k+1)) / D(k);
end
end

内容的提问来源于stack exchange，提问作者Web Inar