波浪模拟过程中大规模计算的优化与加速

有一个任务：模拟海面上的波浪（3D 情况）。我遇到的问题是，当模拟大片海洋（1 公里 x 1 公里）的波浪时，程序需要很长时间才能执行。事实上它根本没有执行。因为在等待结果的半个小时里，它从未被计算过。

我的代码示例：

clear,clc
%% Data 
x = 0:1:100; % coordinates х [m]
y = 0:1:100; % coordinates y [m]
g = 9.81; % gravitational constant [m/s^2]
speed = 5; % wind velocity [m/s]
w0 = (g / speed); % norm.frequency [Hz]
dw = 0.1; % frequency step [rad/s]
w = 0.8:dw:11.1; % frequency [rad/s] 
dtt = pi / 18; % angular step [rad]
theta = 0:dtt:pi; % direction angles, angles between the wavevector & coordintae axis [rad]

%% P-M spectrum, Frequency-Angular spectrum & Amplitude
Psi = 8.1e-3 .* ((w/w0).^(-5)) .* exp((-0.74) ./ ((w/w0).^(4))); % P-M spectrum [none]
Phi = ((speed)^(5)/g^(3)) * Psi; % self-similar spectrum [s*m^2]
Sw = Phi / 2; % frequency spectrum [s*m^2]

St = cos(theta).^(4); % angular spectrum [none]
Norm = trapz(dtt, St); % norm.coefficient [none]
Swt = Sw .* St'; % frequency-angular spectrum [s*m^2]

eta0 = sqrt((Swt * dw * dtt) ./ Norm); % amplitude [m]

figure(1);
subplot(2,1,1)
plot(w, Psi);
title('$$\Psi$$($$\omega$$) - P-M spectrum', 'Interpreter', 'LaTex');
xlabel('\omega [rad/s]');
ylabel('\Psi [none]');
grid on;
subplot(2,1,2)
plot(w, Swt); 
title('$$S(\omega , \theta)$$($$\omega$$) - frequency-angular spectrum', 'Interpreter', 'LaTex');
xlabel('\omega [rad/s]');
ylabel('S(\omega,\theta) [s*m^2]');
grid on;

%% Setting the initial phase parameter
phase = 2*pi*rand(length(theta),length(w)); %% random initial phase ranging from 0 to 2pi [rad]

%% Surface Waving [Linear, 3D (eta & x,y)] at different harmonics & random phase (at one moment in time), different directions of the wavevector (multiple angles)
t = 0; % time moment [s]
Kabs = (w.^2) / g; % wavevector module [rad/m]
Kx = Kabs .* cos(theta)'; % projection of the wavevector onto the x-axis [rad/m]
Ky = Kabs .* sin(theta)'; % projection of the wavevector onto the y-axis [rad/m]

eta = zeros(length(x),length(y),length(theta),length(w)); % reserving space for calculation results
tic
for i = 1:length(x)
    for j = 1:length(y)
        eta(i,j,:,:) = eta0 .* cos(w * t - Kx .* i - Ky .* j + phase);
    end
end
toc
% sum(sum(eta,4),3) - double sum of eta over all harmonics (frequencies) and wavevector directions (angles theta),
% where '4' и '3' summation indicator for variable frequency and angle
etaW = sum(eta,4);
etaWA = sum(etaW,3);

figure(2)
surf(x,y,etaWA);
title('\eta(x,y) - surface waving');
xlabel('x [m]');
ylabel('y [m]');
zlabel('\eta [m]');
cbar = colorbar;
cbar.Label.String = '\eta [m]';
grid on
shading flat

我能够使用的代码优化方法之一是创建一个“空”的 4D 数组（零的 4D 数组）eta = zeros(length(x),length(y),length(theta),length(w));，在执行循环后，计算结果将填充到该数组中：

eta = zeros(length(x),length(y),length(theta),length(w)); % reserving space for calculation results
tic
for i = 1:length(x)
    for j = 1:length(y)
        eta(i,j,:,:) = eta0 .* cos(w * t - Kx .* i - Ky .* j + phase);
    end
end
toc

然后我按频率和角度变量总结结果：

etaW = sum(eta,4);
etaWA = sum(etaW,3);

这样就提前为结果准备好了地方。这有点帮助。例如，x = 0:1:100; y = 0:1:100;使用此方法对一个区域（100 mx 100 m）的代码执行时间为 0.7 秒（不使用则为 3.9 秒）。对于一个区域x = 0:1:500; y = 0:1:500;（500 mx 500 m），执行时间约为 19 秒（不使用则...我不知道，我没有等待代码执行，但结果发现它很长）。但是，对于一个区域x = 0:1:1000; y = 0:1:1000;（1000 mx 1000 m），我很长时间都没有得到想要的结果（感觉根本得不到）。

在我的情况下，还有其他方法可以实现所需的结果并优化我的代码，以便它可以应对如此规模的数据（同时，无需改变数组中的步骤）？

注意：我的电脑有 16 GB 的 RAM。我的第二台电脑（我最初用来进行计算的那台电脑）在执行程序时完全挂断了，所以我不得不“手动”重启它。所以我猜它的 RAM 更少了。