индексы нижнего индекса должны быть либо положительными целыми числами меньше 2 ^ 31, либо логическими - PullRequest
0 голосов
/ 08 мая 2020

SOS У меня постоянно возникают ошибки при решении l oop методом конечных разностей.

Я либо получаю следующую ошибку, когда начинаю с 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' );

1 Ответ

1 голос
/ 09 мая 2020

Взгляните на отредактированный код, который я очистил для вас выше, и изучите его. Не недооценивайте важность чистого, читаемого кода при поиске ошибок. Это сэкономит вам больше времени, чем будет стоить. Особенно через неделю, когда вам нужно будет вернуться к этому коду, и вы совсем не будете вспоминать, что пытались сделать.

Теперь о ваших ошибках. (все ссылки на строки относятся к очищенному коду выше)

Сценарий 1:

В строке 29 вы инициализируете u как одно значение.

Если вы начнете свой l oop в строке 35, начиная с i = 2, то, как только вы попытаетесь выполнить u(i, j+1), т.е. u(2,2) в следующей строке, октава будет жаловаться, что вы пытаетесь для индексации второй строки в массиве, который пока содержит только одну строку. (на самом деле, то же самое применимо и к j на этом этапе, поскольку на этом этапе у вас также есть только один столбец)

Сценарий 2:

Я предполагаю, что Второй сценарий был опечаткой, и вы хотели сказать i = 1 : N. Если вы начинаете с i=1 в l oop, посмотрите на строку 38: вы пытаетесь получить элемент u(i-1, j), т.е. u(0,1). Поэтому октава будет жаловаться, что вы пытаетесь получить элемент ноль , но в массивах октав начинается с единица , а ноль не определен. Попытка получить доступ к любому массиву с нулем приведет к ошибке, которую вы видите (попробуйте в терминале!).

UPDATE

Кроме того, теперь, когда код clean, вы можете обнаружить еще одну ошибку, о которой октава предупредит вас, если вы попытаетесь запустить код.

Посмотрите на строку 26. В отрезке elseif НЕТ условия, поэтому октава ищет следующую заявление в качестве условия проверки.

Это означает, что условие elseif всегда будет успешным, пока результат u0 = 1-x не равен нулю.

Это явно ошибка. Либо вы забыли указать условие для elseif, либо, что более вероятно, вы просто хотели сказать else, а не elseif.

...