Индекс превышает ошибку размеров матрицы в методе Рунге-Кутты: Matlab

Question

Я пытаюсь сделать временный шаговый код с использованием метода Runge-Kutta 4-го порядка, но я столкнулся с проблемами, которые правильно индексируют одно из моих значений. Мой код:

clc; 
clear all; 
L = 32; M = 32; N = 32;      % No. of elements 
Lx = 2; Ly = 2; Lz = 2;      % Size of each element 
dx = Lx/L; dy = Ly/M; dz = Lz/N;   % Step size 
Tt = 1; 
t0 = 0;          % Initial condition 
T = 50;          % Final time 
dt = (Tt-t0)/T;        % Determining time step interval 

% Wave characteristics 
H = 2;          % Wave height 
a = H/2;         % Amplitude 
Te = 6;          % Period 
omega = 2*pi/Te;       % Wave rotational frequency 
d = 25;          % Water depth 
x = 0;          % Location of cylinder axis 

u0(1:L,1:M,1:N,1) = 0;      % Setting up solution space matrix (u values) 
v0(1:L,1:M,1:N,1) = 0;      % Setting up solution space matrix (v values) 
w0(1:L,1:M,1:N,1) = 0;      % Setting up solution space matrix (w values) 

[k,L] = disp(d,omega);      % Solving for k and wavelength using Newton-Raphson function 
%u = zeros(1,50); 
%v = zeros(1,50); 
%w = zeros(1,50); 
time = 1:1:50; 

for t = 1:T 
    for i = 1:L 
     for j = 1:M 
      for k = 1:N 
       eta(i,j,k,t) = a*cos(omega*time(1,t); 
       u(i,j,k,1) = u0(i,j,k,1);  
       v(i,j,k,1) = v0(i,j,k,1); 
       w(i,j,k,1) = w0(i,j,k,1); 
       umag(i,j,k,t) = a*omega*(cosh(k*(d+eta(i,j,k,t))))/sinh(k*d); 
       vmag(i,j,k,t) = 0; 
       wmag(i,j,k,t) = -a*omega*(sinh(k*(d+eta(i,j,k,t))))/sinh(k*d); 
       uRHS(i,j,k,t) = umag(i,j,k,t)*cos(k*x-omega*t); 
       vRHS(i,j,k,t) = vmag(i,j,k,t)*sin(k*x-omega*t); 
       wRHS(i,j,k,t) = wmag(i,j,k,t)*sin(k*x-omega*t); 

       k1x(i,j,k,t) = dt*uRHS(i,j,k,t); 
       k2x(i,j,k,t) = dt*(0.5*k1x(i,j,k,t) + dt*uRHS(i,j,k,t)); 
       k3x(i,j,k,t) = dt*(0.5*k2x(i,j,k,t) + dt*uRHS(i,j,k,t)); 
       k4x(i,j,k,t) = dt*(k3x(i,j,k,t) + dt*uRHS(i,j,k,t)); 
       u(i,j,k,t+1) = u(i,j,k,t) + (1/6)*(k1x(i,j,k,t) + 2*k2x(i,j,k,t) + 2*k3x(i,j,k,t) + k4x(i,j,k,t)); 
       k1y(i,j,k,t) = dt*vRHS(i,j,k,t); 
       k2y(i,j,k,t) = dt*(0.5*k1y(i,j,k,t) + dt*vRHS(i,j,k,t)); 
       k3y(i,j,k,t) = dt*(0.5*k2y(i,j,k,t) + dt*vRHS(i,j,k,t)); 
       k4y(i,j,k,t) = dt*(k3y(i,j,k,t) + dt*vRHS(i,j,k,t)); 
       v(i,j,k,t+1) = v(i,j,k,t) + (1/6)*(k1y(i,j,k,t) + 2*k2y(i,j,k,t) + 2*k3y(i,j,k,t) + k4y(i,j,k,t)); 
       k1z(i,j,k,t) = dt*wRHS(i,j,k,t); 
       k2z(i,j,k,t) = dt*(0.5*k1z(i,j,k,t) + dt*wRHS(i,j,k,t)); 
       k3z(i,j,k,t) = dt*(0.5*k2z(i,j,k,t) + dt*wRHS(i,j,k,t)); 
       k4z(i,j,k,t) = dt*(k3z(i,j,k,t) + dt*wRHS(i,j,k,t)); 
       w(i,j,k,t+1) = w(i,j,k,t) + (1/6)*(k1z(i,j,k,t) + 2*k2z(i,j,k,t) + 2*k3z(i,j,k,t) + k4z(i,j,k,t)); 
       a(i,j,k,t+1) = ((u(i,j,k,t+1))^2 + (v(i,j,k,t+1))^2 + (w(i,j,k,t+1))^2)^0.5; 
      end 
     end 
    end 
end 

На данный момент значения, кажется, будет хорошо для первой итерации, но тогда у меня есть ошибка Index exceeds matrix dimension в строке расчета eta. Я понимаю, что я неправильно индексирую значение eta, но не уверен, как это исправить.

Моя цель - обновить значение eta для каждого цикла t, а затем использовать это новое значение eta для остальных вычислений.

Я по-прежнему новичок в программировании и стараюсь понять индексацию, особенно в 3 или 4-мерных матрицах, и будет очень признателен за любые советы по правильному вычислению этого значения.

Заранее благодарим за любой совет!

Answer 1

Вы объявляете

time = 1:1:50; 

, который только вектор-строка, но доступ к нему здесь

eta(i,j,k,t) = a*cos(omega*time(i,j,k,t)); 

, как если бы это был массив с 4-х измерениях.

Для корректного доступа к элементу x из time вам нужно использовать синтаксис

time(1,x); 

(как это массив 1 х 50)