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

Question

Я пытаюсь сделать код с пошаговым шагом времени, используя метод Рунге-Кутты 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)