matlab写涡度方程各项的代码问题

小蛇

                               
登录/注册后可看大图

因为需要涡度方程，以家园里的分享grads代码想写一个matlab版本的，花了一天时间写了一个M项，但调试了半天不知为什么好几项都是空值，求出的V2项也仅仅是一个值而不是设想的一个数组，报错和代码如下，求大神解答。。。本科生跪谢。。。。。

% 读取数据
nc_path = '202007.nc';  % 数据路径
ncid = netcdf.open(nc_path, 'nowrite');
lon = netcdf.getVar(ncid, netcdf.inqVarID(ncid, 'longitude'));
lat = netcdf.getVar(ncid, netcdf.inqVarID(ncid, 'latitude'));
level = netcdf.getVar(ncid, netcdf.inqVarID(ncid, 'level'));
time = netcdf.getVar(ncid, netcdf.inqVarID(ncid, 'time'));
u = netcdf.getVar(ncid, netcdf.inqVarID(ncid, 'u'));
v = netcdf.getVar(ncid, netcdf.inqVarID(ncid, 'v'));
vo = netcdf.getVar(ncid, netcdf.inqVarID(ncid, 'vo'));
netcdf.close(ncid);
% 计算涡度方程各项
[lon1, lat1]=meshgrid(lon,lat);
for x = 1:length(lon1);  % 经度范围
     for y = 1:length(lat1);  % 纬度范围
        z0 = 1;  % 垂直层数
        t0 = 1;  % 时间步长
        fi = lat(y) * pi / 180;  % 纬度转弧度
        beta = 2 * cos(fi) * 7.29e-5;  % beta参数
        % 计算dvo_dx和dvo_dy
        % 计算dvo_dx和dvo_dy
        dx = diff(lon(x) * pi / 180);  % 计算经度间隔
        dy = diff(lat(y) * pi / 180);  % 计算纬度间隔
        dvo_dx = double(diff(vo(x, y, z0, t0), 1, 1)) ./ double(dx) * 6371000 * cos(fi);  % x方向上的vo梯度
        dvo_dy = double(diff(vo(x, y, z0, t0), 1, 2)) ./ double(dy) * 6371000;  % y方向上的vo梯度
        
        V1 = -double(u(x, y, z0, t0)) .* dvo_dx - double(v(x, y, z0, t0)) .* dvo_dy;
        V4 = double(v(x, y, z0, t0)) .* beta;
        V1_all(x, y) = V1;  % 将计算结果存储到数组中
        V4_all(x, y) = V4;
    end
end

                               
登录/注册后可看大图

                               
登录/注册后可看大图

Sigmas

1. % 读取数据
   
2. nc_path = '202007.nc';  % 数据路径
   
3. ncid = netcdf.open(nc_path, 'nowrite');
   
4. lon = netcdf.getVar(ncid, netcdf.inqVarID(ncid, 'longitude'));
   
5. lat = netcdf.getVar(ncid, netcdf.inqVarID(ncid, 'latitude'));
   
6. level = netcdf.getVar(ncid, netcdf.inqVarID(ncid, 'level'));
   
7. time = netcdf.getVar(ncid, netcdf.inqVarID(ncid, 'time'));
   
8. u = netcdf.getVar(ncid, netcdf.inqVarID(ncid, 'u'));
   
9. v = netcdf.getVar(ncid, netcdf.inqVarID(ncid, 'v'));
   
10. vo = netcdf.getVar(ncid, netcdf.inqVarID(ncid, 'vo'));
    
11. netcdf.close(ncid);
    
12. 
    
13. % 计算涡度方程各项
    
14. [lon1, lat1] = meshgrid(lon, lat);
    
15. V1_all = zeros(length(lon), length(lat));
    
16. V4_all = zeros(length(lon), length(lat));
    
17. 
    
18. for x = 1:length(lon1);  % 经度范围
    
19.     for y = 1:length(lat1);  % 纬度范围
    
20.         z0 = 1;  % 垂直层数
    
21.         t0 = 1;  % 时间步长
    
22.         fi = lat(y) * pi / 180;  % 纬度转弧度
    
23.         beta = 2 * cos(fi) * 7.29e-5;  % beta参数
    
24.         
    
25.         % 计算dvo_dx和dvo_dy
    
26.         dx = (lon(2) - lon(1)) * pi / 180;  % 计算经度间隔
    
27.         dy = (lat(2) - lat(1)) * pi / 180;  % 计算纬度间隔
    
28.         dvo_dx = double(diff(vo(x, :, z0, t0), 1, 2)) ./ double(dx) * 6371000 * cos(fi);  % x方向上的vo梯度
    
29.         dvo_dy = double(diff(vo(:, y, z0, t0), 1, 1)) ./ double(dy) * 6371000;  % y方向上的vo梯度
    
30.         % 补充缺失的最后一个元素
    
31.         dvo_dx = [dvo_dx, dvo_dx(end)];
    
32.         dvo_dy = [dvo_dy; dvo_dy(end)];
    
33. 
    
34.         V1 = -double(u(x, y, z0, t0)) .* dvo_dx - double(v(x, y, z0, t0)) .* dvo_dy;
    
35.         V4 = double(v(x, y, z0, t0)) .* beta;
    
36.         V1_all(x, y) = V1;  % 将计算结果存储到数组中
    
37.         V4_all(x, y) = V4;
    
38.     end
    
39. end

复制代码

盲写的，我也不知道能不能解决你的问题，如果有报错就更好了

小蛇

Sigmas 发表于 2023-3-27 00:51
盲写的，我也不知道能不能解决你的问题，如果有报错就更好了

非常感谢您的回答，我按照您的修改了一下，有v1_all数组了，但依然是空值，不会附图，报错如下：
错误使用  ./
对于此运算，数组的大小不兼容。

出错 V1ANDV4 (第 25 行)
        dvo_dx = double(diff(vo(x, :, z0, t0), 1, 1)) ./ double(dx) * 6371000 * cos(fi);  % x方向上的vo梯度

wjy_ecnu

for i=1:length(lon1)
for j=1:length(lat1)
这里是有问题的
应该是[m,n]=size(lon1)
for j=1:m
for j=1:n
xxxxxx
end
end   
这样就能避免你说的空值得问题

wjy_ecnu

另外 实际上matlab应该不需要这么多层for循环   式子的每一项用一行应该就能搞定 只是那个计算相邻格点距离的时候 可能用循环更好理解一点  用向量化的方法计算格点距离 可能会看起来不好理解

小蛇

wjy_ecnu 发表于 2023-3-27 10:09
另外 实际上matlab应该不需要这么多层for循环   式子的每一项用一行应该就能搞定 只是那个计算相邻格点距离 ...

不用for循环怎么求那么多格点上的值呢？

Sigmas

小蛇 发表于 2023-3-27 10:05
非常感谢您的回答，我按照您的修改了一下，有v1_all数组了，但依然是空值，不会附图，报错如下：
错误使 ...

1. dvo_dx = double(diff(vo(x, :, z0, t0), 1, 2)) * 6371000 * cos(fi) ./ double(dx);  % x方向上的vo梯度
   
2. dvo_dy = double(diff(vo(:, y, z0, t0), 1, 1)) * 6371000 ./ double(dy);  % y方向上的vo梯度

复制代码

改一下计算顺序试试呢？
另外以我浅薄的理解也许可以稍微回答一下你对于下面老哥“如何不用for循环”的疑惑，在这个实例中用for循环确实更直观，但改用其他方式运算也许速度可以更快并且更简洁（比如用wjy老哥说到的矩阵向量化）

wjy_ecnu

散度 divergence  涡度 curl  梯度gradient  梯度注意matlab会自动改变x和y的顺序   但是这些函数不直观
如果研究一下 仍然没有理解 还是用for循环吧  毕竟循环逻辑清楚

小蛇

wjy_ecnu 发表于 2023-3-27 15:04
散度 divergence  涡度 curl  梯度gradient  梯度注意matlab会自动改变x和y的顺序   但是这些函数不直观
...

好的，谢谢您的建议。

小蛇

Sigmas 发表于 2023-3-27 12:36
改一下计算顺序试试呢？
另外以我浅薄的理解也许可以稍微回答一下你对于下面老哥“如何不用for循环” ...

我再试试，非常感谢