user3460758
user3460758

Reputation: 987

Index exceeds matrix dimensions error in Runge-Kutta method: Matlab

I'm trying to make a time stepping code using the 4th order Runge-Kutta method but am running into issues indexing one of my values properly. My code is:

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

At the moment, the values seem to be fine for the first iteration but then I have the error Index exceeds matrix dimension in the line calculating eta. I understand that I am not correctly indexing the eta value but am not sure how to correct this.

My goal is to update the value of eta for each loop of t and then use that new eta value for the rest of the calculations.

I'm still quite new to programming and am trying to understand indexing, especially in 3 or 4 dimensional matrices and would really appreciate any advice in correctly calculating this value.

Thanks in advance for any advice!

Upvotes: 0

Views: 148

Answers (1)

mathematician1975
mathematician1975

Reputation: 21351

You declare

time = 1:1:50;

which is just a row vector but access it here

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

as if it were an array with 4 dimensions.

To correctly access element x of time you need to use syntax

time(1,x);

(as it is a 1 x 50 array)

Upvotes: 1

Related Questions