векторизовать 2Д градиент с пространственной вариацией закромах


Следующий код принимает некоторые значения ssh и решает уравнения геострофического движения (слайд 8 здесь).

Основная часть кода является вычисление частных производных ssh. В частности, дискретные различия должны быть умножены на коэффициент, который зависит от y. Это можно легко сделать с линейный код:

close all
clearvars
clc

%define grid
x=linspace(120,145,20);
y=linspace(30,45,20);
[x, y]=meshgrid(x,y);
x = x'; 
y = y';

%define constants
R = 6371000; % Earth radius

g=9.806-.5*(9.832-9.780)*cos(2*y*pi/180); % gravity

omega = 2*pi/(24*60*60); % Earth rotation angle velocity [s]
f = 2*omega*sind(y); %Coriolis force coefficients


%data
ssh = (exp(-((x-130).^2/20)).*(exp(-(y-35).^2/7)))*1e6; % Sea surface height in each point


%Calculate geostrophic current
u=zeros(size(ssh));
v=zeros(size(ssh)); 
for i=2:size(x,1)-1
    for j=2:size(y,2)-1
        dx(i,j) = (x(i+1,j)-x(i-1,j)) *(R*cosd(y(i,j))*pi/180);
        dy(i,j) = (y(i,j+1)-y(i,j-1)) *(R*pi/180);

        u(i,j) = -g(i,j)/f(i,j) *(ssh(i,j+1)-ssh(i,j-1)) /dy(i,j);
        v(i,j) =  g(i,j)/f(i,j) *(ssh(i+1,j)-ssh(i-1,j)) /dx(i,j);
    end
end

figure
pcolor(x,y,ssh)
shading flat
hold on
quiver(x,y,u,v,2,'k')
title('Geostrophic current [m/s]','fontweight','bold')
xlabel('longitude','fontweight','bold')
ylabel('latitude','fontweight','bold')
set(gcf,'color','w')

Выход: first code

Однако, у меня возникли проблемы, чтобы векторизовать код.

Я пытался использовать gradient функция следующим образом:

%%%%===== vectorized code =====%%%%

dx2 = x  .*(R*cosd(y)*pi/180); %x-position matrix
dy2 = y  *(R*pi/180); %y-position matrix

[dsshdy,dsshdx] = gradient(ssh, dy2,dx2);

u2 = -g./f .*dsshdy; 
v2 =  g./f .*dsshdx; 

figure;hold on
pcolor(x,y, ssh)
shading flat
hold on
quiver(x,y, u2,v2, 2,'k')
title('Geostrophic current 2','fontweight','bold')
xlabel('longitude','fontweight','bold')
ylabel('latitude','fontweight','bold')
set(gcf,'color','w')

Выход: second code

Однако, это не так (я думаю), так как функция градиента не принимает в качестве входных данных матрицы значений интервала. Как следствие, код как-то вычисляет дифференциалы, которые являются слишком большой, и стрелки не visibile.

Как я могу векторизовать такую проблему без повторного введения цикла с учетом изменения dx С y?



114
1
задан 6 февраля 2018 в 09:02 Источник Поделиться
Комментарии
1 ответ

Ваш f и ssh уже векторизованных, вы можете сделать то же самое довольно тривиально с u и v также. Нет ничего хитрого происходит в цикле. Процесс просто удалить forоставляя задание i=2:size(x,1)-1. И заменить все матрицы умножение и деление на поэлементного умножения и деления (.* и ./). Эта листья:

%Calculate geostrophic current
u = zeros(size(ssh));
v = zeros(size(ssh));
i = 2:size(x,1)-1;
j = 2:size(y,2)-1;
dx(i,j) = (x(i+1,j)-x(i-1,j)) .* (R*cosd(y(i,j))*pi/180);
dy(i,j) = (y(i,j+1)-y(i,j-1)) .* (R*pi/180);
u(i,j) = -g(i,j) ./ f(i,j) .* (ssh(i,j+1)-ssh(i,j-1)) ./ dy(i,j);
v(i,j) = g(i,j) ./ f(i,j) .* (ssh(i+1,j)-ssh(i-1,j)) ./ dx(i,j);

Вы можете затем сделать небольшое упрощение, dx и dy не нуждаются в индексации, так как вы используете ту же часть, что вам назначить:

dx = (x(i+1,j)-x(i-1,j)) .* (R*cosd(y(i,j))*pi/180);
dy = (y(i,j+1)-y(i,j-1)) .* (R*pi/180);
u(i,j) = -g(i,j) ./ f(i,j) .* (ssh(i,j+1)-ssh(i,j-1)) ./ dy;
v(i,j) = g(i,j) ./ f(i,j) .* (ssh(i+1,j)-ssh(i-1,j)) ./ dx;

1
ответ дан 9 февраля 2018 в 06:02 Источник Поделиться