0

我正在尝试使用 4 阶 Runge-Kutta 方法制作时间步进代码，但我遇到了正确索引我的值之一的问题。我的代码是：

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 值，但不确定如何更正此问题。

我的目标是为每个 t 循环更新 eta 的值，然后将新的 eta 值用于其余计算。

我对编程还是很陌生，并且正在尝试理解索引，尤其是在 3 或 4 维矩阵中，并且非常感谢有关正确计算此值的任何建议。

提前感谢您的任何建议！

4

1 回答 1

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 x 50 阵列）

于 2015-05-01T14:05:51.003 回答