Запись векторной суммы в MATLAB - PullRequest
Новая
       30

Запись векторной суммы в MATLAB

0 голосов
/ 11 марта 2020

Предположим, у меня есть функция phi(x1,x2)=k1*x1+k2*x2, которую я оценил по сетке, где сетка представляет собой квадрат с границами -100 и 100 по осям x1 и x2 с некоторым размером шага, скажем h=0.1. Теперь я хочу вычислить эту сумму по сетке, с которой я борюсь:

enter image description here

Что я пытался: 

clear all
close all
clc
D=1; h=0.1;
D1 = -100;
D2 = 100;
X = D1 : h : D2;
Y = D1 : h : D2;
[x1, x2] = meshgrid(X, Y);
k1=2;k2=2;
phi = k1.*x1 + k2.*x2;
figure(1)
surf(X,Y,phi)
m1=-500:500;
m2=-500:500;
[M1,M2,X1,X2]=ndgrid(m1,m2,X,Y)
sys=@(m1,m2,X,Y) (k1*h*m1+k2*h*m2).*exp((-([X Y]-h*[m1 m2]).^2)./(h^2*D))
sum1=sum(sys(M1,M2,X1,X2))

Matlab говорит об ошибке в ndgrid, есть идеи, как мне это закодировать?

MATLAB показывает: 

Error using repmat
Requested 10001x1001x2001x2001 (298649.5GB) array exceeds maximum array size preference. Creation of arrays greater
than this limit may take a long time and cause MATLAB to become unresponsive. See array size limit or preference
panel for more information.

Error in ndgrid (line 72)
        varargout{i} = repmat(x,s);

Error in new_try1 (line 16)
[M1,M2,X1,X2]=ndgrid(m1,m2,X,Y)

1 Ответ

2 голосов
/ 12 марта 2020

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

Чтобы получить значение M (x1, x2) для некоторого заданного значения ( x1, x2), вы должны вычислить эту сумму по Z2. Конечно, используя числовой набор инструментов, такой как MATLAB, вы можете рассчитывать только на некоторый конечный диапазон Z2. В этом случае, поскольку (x1, x2) охватывает диапазон [-100,100] x [-100,100], а h = 0,1, отсюда следует, что mh охватывает диапазон [-1000, 1000] x [-1000, 1000]. Пример: m = (-1000, -1000) дает вам mh = (-100, -100), который является левым нижним углом вашего домена. Так что действительно, фи (мч) - это просто фи (х1, х2), оцениваемая по всем вашим дискретным точкам.

Кроме того, поскольку вам необходимо вычислить |x-hm|^2, вы можете рассматривать x = x1 + i x2 как комплексное число, чтобы использовать функцию abs MATLAB. Если бы вы строго работали с векторами, вам пришлось бы использовать norm, что тоже нормально, но более многословно. Таким образом, для некоторого заданного x=(x10, x20) вы бы вычислили x-hm по всей дискретной плоскости как (x10 - x1) + i (x20 - x2).

Наконец, вы можете вычислить 1 член M за раз: 

D=1; h=0.1;
D1 = -100;
D2 = 100;
X = (D1 : h : D2); % X is in rows (dim 2)
Y = (D1 : h : D2)'; % Y is in columns (dim 1)
k1=2;k2=2;
phi = k1*X + k2*Y;

M = zeros(length(Y), length(X));

for j = 1:length(X)
    for i = 1:length(Y)
        % treat (x - hm) as a complex number
        x_hm = (X(j)-X) + 1i*(Y(i)-Y); % this computes x-hm for all m
        M(i,j) = 1/(pi*D) * sum(sum(phi .* exp(-abs(x_hm).^2/(h^2*D)), 1), 2);
    end
end

Кстати, это вычисление занимает довольно много времени. Вы можете рассмотреть возможность увеличения h, уменьшения D1 и D2 или изменения всех трех из них.

...