所以我做了这个练习:
利用公式(18),实现Gauss-Legendre方法。附加函数 [t, w] = gauleg(n) 返回给定多项式次数 n 的节点值 ti 和权重 wi。
这是我的代码:
% Function to generate nodes and weights for Gauss-Legendre quadrature
function [t, w] = gauleg(n)
% Create a matrix related to Legendre polynomials
beta = 0.5 ./ sqrt(1 - (2 * (1:n-1)).^(-2)); % Calculate parameters
T = diag(beta, 1) + diag(beta, -1); % Tridiagonal matrix
% Compute eigenvalues and eigenvectors
[V, D] = eig(T); % Eigenvalue decomposition
% Nodes are the eigenvalues
t = diag(D); % Nodes
[t, idx] = sort(t); % Sort nodes
% Weights are the square of the first element of each eigenvector, multiplied by 2
w = 2 * (V(1, idx)).^2; % Compute weights based on sorting index
end
% Function to compute integrals using Gauss-Legendre quadrature
function I = f_gaussLegendre(f, a, b, n)
% Get the nodes and weights
[t, w] = gauleg(n);
% Midpoint and half of the interval's length
c = (a + b) / 2; % Midpoint of the interval
d = (b - a) / 2; % Half of the interval's length
% Compute function values at the correct points
F = arrayfun(@(ti) f(c + d * ti), t); % Compute function values at the nodes
weighted_sum = sum(w .* F); % Sum the weighted function values
% Compute the final integral result
I = d * weighted_sum; % Final integration result
end
% Test code
% Integration interval
a = 0;
b = 2;
% Function to integrate
f = @(x) x.^3 + 2*x.^2 + 1;
% Compute integral using Gauss-Legendre method
n_values = [2, 4, 6]; % Example values for n
results = zeros(size(n_values)); % Allocate space for results
% Loop through different values of n
for i = 1:length(n_values)
results(i) = f_gaussLegendre(f, a, b, n_values(i));
end
% Display the results
disp("Results from Gauss-Legendre method:");
disp("n | result");
disp([n_values; results]'); % Display the values for n and their corresponding results
这段代码的输出是
error: =: nonconformant arguments (op1 is 1x1, op2 is 1x2)
Results from Gauss-Legendre method:
n | result
2 0
4 0
6 0
“我不知道错误在哪里,也不知道如何修复它,因为 Jupyter Notebook 不显示哪一行有错误。
首先检查单个计算总是一个好主意,而不是立即通过循环进行整个计算。跑步
f_gaussLegendre(f, a, b, n_values(1))
表明您的函数返回一个向量
[11.333 11.333]
,它不适合标量变量 result(1)
,因此 Octave 会报告上述错误。
%debug
魔法命令是否适用于 Octave 的 Jupyter,但如果不能,您始终可以恢复到良好可信的古老方式,只需打印出值进行调试。通过 打印出
w
、
F
和
weighed_sum
变量的大小
...
weighted_sum = sum(w .* F); % Sum the weighted function values
disp(size(w)); disp(size(F)); disp(size(weighted_sum));
...
显示尺寸分别为
[1 2]
、[2 1]
和 [1 2]
。很明显,存在问题,因为乘法 w .* F
会产生 2x2 矩阵,而 sum
将简单地对其列求和。也许你只是想在那里做 w*F
,或者可能某个地方缺少转置。