Varga's questions

有一组 excel 文件（10 个文件），每个文件都有一个表格：2 列（参数名称、参数值）和 12 行（参数数量及其值）。需要计算单个参数的 10 个变量之间的平均值。下面有一个程序可以实现此目的。但是，我注意到他执行计算时的误差为万分之一，他低估了它。例如，结果如下：

clc, clear;
 
% Folder with Excel files
folder = 'The path to the folder with 10 excel files'; % the path to the folder
files = dir(fullfile(folder, '*.xlsx')); % we get a list of all excel files
 
% Initialize the vector to store the average values
average_values = [];
 
% Loop through all files
for i = 1:length(files)
    % The full path to the file
    file_path = fullfile(folder, files(i).name);
 
    % Reading data from an Excel file
    opts = detectImportOptions(file_path, 'VariableNamingRule', 'preserve'); 
    data = readtable(file_path, opts);
 
    % Check if there is at least one numeric column
    if any(varfun(@isnumeric, data, 'OutputFormat', 'uniform'))
        % We extract only numeric data
        numeric_data = data{:, varfun(@isnumeric, data, 'OutputFormat', 'uniform')}; 
 
        % We calculate the average values for the rows, ignoring NaN
        avg_values_per_row = mean(numeric_data, 2, 'omitnan');
 
        % Saving the average values in an array
        average_values = avg_values_per_row;
    else
        warning('There are no numeric columns in the "%s" file.', files(i).name);
    end
end
 
% Displaying average values
disp('Average values by line:');
disp(average_values);
 
% Displaying the number of rows
disp('Number of lines:');
disp(length(average_values));

如您所见，对于参数编号 1、4、6、8，低估值从 0.0001 到 0.0007 个单位（对于第 8 个参数，低估值已经是该数字的千分之一）。最有趣的是，如果您手动将每个参数的 10 个值输入meanMatlab 中的函数并进行计算，它将给出正确的结果，而不会出现这种低估（使用第 8 个参数的示例）：

>> mean([0.952380002    0.980380905 0.97087088  0.95237715  0.952379051 0.970871364 0.980383351 0.925921306 0.961536062 0.952370972
])

ans =

    0.9599

代码中哪里可能存在问题以及如何修复它以便给出准确的结果？（附件包含一个带有示例文件的文件夹： https: //drive.google.com/file/d/1COPpBq2mU8Ww62DRvbOKTtr_aC-LsCIL/view? usp=sharing ）

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

我有 101 个 .csv 文件，每个文件有 10001 对值（图中此类文件的一个示例是频谱功率密度对频率的依赖性 -示例）。我需要以某种方式对每个文件中第二列的值进行平方，然后对这些值进行平均（平均发生在一行号的值之间。即，我们从所有 101 个文件中获取一行第二列的平方值，对它们进行平均，对于所有其他行也是如此。换句话说，依赖性应该在频谱功率密度上进行平均。

如何才能做到这一点？

到目前为止，目前我只能处理一个文件，同时，一个文件的依赖关系图（没有平方）看起来就像一场噩梦。

import pandas as pd
import matplotlib.pyplot as plt

data = pd.read_csv(r'C:\Users\Den\Desktop\SPD\tek0000MTH.csv', sep = ';')


x = data['W'].to_numpy()
y = data['S'].to_numpy()

plt.xlabel('w') 
plt.ylabel('S') 
plt.title('SPD') 
plt.plot(x, y)
plt.grid() 
plt.show()

阴谋

同时，我无法想象如何执行上述所有任务，例如平方和平均，同时处理 101 个文件。