时间戳

首页 > 解决方案 > 下标索引必须是小于 2^31 的正整数或逻辑

algorithm - 下标索引必须是小于 2^31 的正整数或逻辑

问题描述

SOS 我在有限差分法的循环求解中不断出错。

当我开始时，我会收到以下错误i = 2 : N：

diffusion: A(I,J): row index out of bounds; value 2 out of bound 1
error: called from
    diffusion at line 37 column 10  % note line change due to edit!

或者，当我这样做时出现以下错误i = 2 : N：

subscript indices must be either positive integers less than 2^31 or logicals
error: called from
    diffusion at line 37 column 10   % note line change due to edit!

请帮忙

  clear all; close all;

% mesh in space
  dx = 0.1;
  x  = 0 : dx : 1;

% mesh in time
  dt = 1 / 50;
  t0 = 0;
  tf = 10;
  t  = t0 : dt : tf;

% diffusivity
  D = 0.5;

% number of nodes
  N = 11;

% number of iterations
  M = 10;

% initial conditions
  if x <= .5 && x >= 0   % note, in octave, you don't need parentheses around the test expression
      u0 = x;
  elseif 
      u0 = 1-x;
  endif

  u = u0;

  alpha = D * dt / (dx^2); 

  for j = 1 : M

      for i = 1 : N
          u(i, j+1) = u(i, j )            ...
                      + alpha             ...
                        * ( u(i-1, j)     ...
                            + u(i+1, j)   ...
                            - 2           ...
                              * u(i, j)   ...
                          )               ;
      end

      u(N+1, j+1) = u(N+1, j)             ...
                    + alpha               ...
                      * (                 ...
                          u(N, j)         ...
                          - 2             ...
                            * u(N+1, j)   ...
                          + u(N, j)       ...
                        )                 ;

    % boundary conditions
      u(0, :) = u0;
      u(1, :) = u1;
      u1      = u0;
      u0      = 0;

  end 


% exact solution with 14 terms

  %k=14   % COMMENTED OUT

  v = (4 / ((k * pi) .^ 2))             ...
      * sin( (k * pi) / 2 )             ...
      * sin( k * pi * x )               ...
      * exp .^ (D * ((k * pi) ^ 2) * t) ;

  exact = symsum( v, k, 1, 14 );
  error = exact - u;

% plot stuff
  plot( t, error );
  xlabel( 'time' );
  ylabel( 'error' );
  legend( 't = 1 / 50' );

标签： algorithmoctave

解决方案

看看我在上面为你清理的编辑代码并研究它。在寻找错误时，不要低估干净、可读的代码的重要性。它将为您节省更多的时间而不是成本。尤其是从现在开始的一周后，当您需要重新访问此代码时，您将完全不记得您正在尝试做什么。

现在关于你的错误。（所有行引用都与上面的清理代码有关）

场景一：

在第 29 行中，您将初始化u为单个值。

如果您从第 35 行开始循环i = 2，那么一旦您尝试执行u(i, j+1)，即u(2,2)在下一行中，octave 就会抱怨您正在尝试索引第二行，该数组目前仅包含一行. （事实上​​，此时 j 同样适用，因为此时您也只有一列）

场景二：

我认为第二种情况是一个错字，你的意思是说i = 1 : N. 如果您从i=1循环开始，请查看第 38 行：您正在尝试获取 element u(i-1, j)，即u(0,1). 因此 octave 会抱怨您试图获取零元素，但在 octave 数组中，数组从1开始并且未定义零。尝试使用零访问任何数组将导致您看到错误（在终端中尝试！）。

更新

此外，既然代码是干净的，您可以发现另一个错误，如果您尝试运行代码，该 octave 会有效地警告您。

看第 26 行。在elseif腿中没有条件，所以 octave 寻找下一条语句作为测试条件。

这意味着elseif只要 u0 = 1-x 的结果不为零，条件将始终成功。

这显然是一个错误。要么您忘记为 设置条件elseif，或者更可能的是，您可能只是想说else，而不是elseif。

推荐阅读