下标指数必须是小于2^31的正整数或逻辑数。

Question

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' );

Answer 1

看一下我上面为你清理的编辑代码，并研究一下.当猎取bug时，不要低估干净、可读的代码的重要性.它将为你节省更多的时间，而不是代价。尤其是一周后，当你需要重新审视这段代码时，你将完全不记得你想做什么。

现在关于你的错误。(所有的行引用都是针对上面的清理后的代码)

场景1:

在第29行，你初始化了 u 作为一个单一的值。

如果你在第35行开始循环，以 i = 2然后，一旦你尝试做 u(i, j+1)，即 u(2,2) 在下一行中，octave会抱怨你在试图将 第二 行，在一个数组中，到目前为止只包含了 一个 行。(事实上，此时j也将同样适用，因为此时你也只有一列)

第二种情况。

我认为第二种情况是一个错字，你的意思是说： i = 1 : N.如果你从 i=1 循环中，然后看看第38行：你试图让元素 u(i-1, j)即： u(0,1). 因此，八度会抱怨你想把 零 元素，但在八度数组中则从 一个 而零没有被定义。试图访问任何带零的数组都会导致你看到的错误（在终端机上试试！）。

更新

另外，现在代码已经干净了，你可以发现另一个错误，如果你尝试运行代码，八度会帮助你警告你。

请看第26行。代码中没有任何条件 elseif 腿，所以八度寻找下一条语句作为测试条件。

这意味着 elseif 只要u0=1-x的结果是非零，条件就会成功。

这显然是一个错误。要么是你忘了将条件放在 elseif或者更有可能的是，你可能只是想说 else而不是 elseif.