有一个任务：模拟海面上的波浪（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 更少了。
