Octave 中的不一致参数

问题描述 投票:0回答:1

所以我做了这个练习:

利用公式(18),实现Gauss-Legendre方法。附加函数 [t, w] = gauleg(n) 返回给定多项式次数 n 的节点值 ti 和权重 wi。

Formula 18

这是我的代码:

% 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 不显示哪一行有错误。

jupyter-notebook jupyter octave integral
1个回答
0
投票

首先检查单个计算总是一个好主意,而不是立即通过循环进行整个计算。跑步

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
,或者可能某个地方缺少转置。

© www.soinside.com 2019 - 2024. All rights reserved.