Построение графика трёхмерной свёртки двух функций в Matlab

@BitMarshal · Регистрация: 15.05.2013

Студворк — интернет-сервис помощи студентам

Доброго времени суток!

Пишу первый раз, поэтому прошу не судить меня строго за плохое оформление.
Задача связана с дипломной работой, поэтому решение необходимо срочно.
А проблема звучит так: необходимо построить на одном рисунке два графика (линий уровня и трёхмерный) трёхмерной свёртки двух функций: функции источника единичной интенсивности f(x,y,z) и функции плотности распределения вещества g(x,y,z).

Функция источника выглядит так:
f(x,y,z) = [1 – cos (arcsin (a/r))] * [1 / (L + R – r)² + 1 / (L + R + r)²] * k * R² / 8, где r = √(x² + y² + z²).
Функция плотности (имеет вид сферы):
g(x,y,z) = z₀ + √(R₀² – (x – x₀)² – (y – y₀)²), где x₀, y₀, z₀ — настраиваемые координаты центра области распределения.
Свёртка этих функций должна иметь такой вид:
I (x,y,z) = ∫∫∫ f(x–x', y–y', z–z') * g(x', y', z') dx' dy' dz'.
В случае моей задачи:
I (x,y,z) = ∫∫∫ [1 – cos(arcsin(a / √((x–x')² + (y–y')² + (z–z')²)))] * [1 / (L + R – √((x–x')² + (y–y')² + (z–z')²))² + 1 / (L + R + √((x–x')² + (y–y')² + (z–z')²))²] * √(R₀² – (x'–x₀)² – (y'–y₀)²) * k * R² / 8 dx' dy' dz'.
Пределы интегрирования — (0, x), (0, y) и (0, z) соответственно.

Проведя численное интегрирование, необходимо построить график поверхности от (x,y,I). Координаты x,y,z должны менятся от -25 до 25 каждая.
В начале программы также должен быть блок вводимых параметров:
a = 0.01 (с клавиатуры), R = 50, L = 50, R₀ = 4, k = 1.
Затем необходимо сделать всё то же самое, но уже в качестве функции плотности будет выступать свёртка двух сфер, т.е.
∫∫ [z₁ + √(R₁² – (x'–x₁)² – (y'–y₁)²)] * [z₂ + √(R₁² – (x–x₂–x')² – (y–y₂–y')²)] dx' dy'.
Проблема заключается в том, что Matlab отказывается проводить интегрирование с помощью команды triplequad (явно интеграл посчитать не удаётся) с переменными пределами.

Прошу мне помочь, так как для меня это очень важно. Matlabэом пользуюсь всего месяц и много чего не знаю. Версия моей программы — R2012b (8.00).

Зосима · 24.05.2013, 17:09

так, хорошо, так что осталось делать? чтот я нить потерял...

@BitMarshal · 24.05.2013, 18:41 **[ТС]**

Кроме самой формулы (это чисто математическая проблема), осталось только понять, как решить проблему суммы слоёв матрицы I(i,j,q).
Для одного слоя z=1 решение задачи я понял. А если надо, чтоб z = 0:n и интегрирование было zz = 0:z плюс суммировать для каждой I(i,j) все значения от z = 0 до z = n, вы предложили обратиться к этому посту. Но вот проблема: там z = const! А мне надо, чтоб суммирование для каждой I(i,j) было по всем z от 0 до n (n = 30).
Как я понимаю, нужен 3-й цикл. Без проблем, но в таком случае я получу 3-мерную матрицу. Вопрос, к сожалению, остаётся: как поэлементно сложить все слои матрицы I(i,j,q) по координате z для каждой точки (i,j) и получить двумерную матрицу Ii(i,j), из которой я уже буду строить график?
Ну, т.е.

Matlab M
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
clc
clear all
% Начальные условия
a = 2;
R = 50; L = 50; k = 1; R0 = 4;
x0 = 15; y0 = 15; z0 = 0;
n = 30; m = 50; l = 10;
% Расчёт
x = 0:n; y = 0:n; z = 0:n;
I = zeros(length(x),length(y),length(z));
for i = 1:length(x)
  for j = 1:length(y)
    for q = 1:length(z)
      [xx yy zz] = meshgrid( linspace(0,x(i),m), linspace(0,y(j),m), linspace(0,z(q),m) );
      r = sqrt( (x(i)-xx).^2 + (y(j)-yy).^2 + (z(q)-zz).^2 );
      fr = z0 + sqrt(R0^2 - (xx-x0).^2 - (yy-y0).^2);
      f = (1 - cos(asin(a./r))).* fr * k * R^2 / (4 * (L + R)^2);
      f( isnan(f)|isinf(f) ) = 0;
      Fdz = trapz(linspace(0,z(q),m), f, 3);
      Fdydz = trapz(yy(:,1,1), Fdz);
      I(i,j,q) = trapz(xx(1,:,1), Fdydz);
    end
  end
end

Что-то я сомневаюсь насчёт правильности 19-й строки.. Ну ладно. Так вот. Я получил матрицу I(i,j,q) размером 31x31x31. Как для каждой точки I(i,j) сложить все значения по слоям?

Зосима · 25.05.2013, 13:41

Ага

кажется дошло

ты правильно рассуждаешь.
Попробуй:
I2 = sum(I, 3);
*здесь 3 - номер размерности, по которой идет суммирование.
Поидее, если я не напортачил, то I2 должно получиться двумерной матричкой. Проверяй

@BitMarshal · 25.05.2013, 14:41 **[ТС]**

Да, ваш способ работает, I(i,j,q) стала I2(i,j). Большое спасибо! График, правда, всё равно отстой, но теперь я точно уверен, что ошибка в формуле.

Построение графика трёхмерной свёртки двух функций в Matlab

БЛАГОДАРЮ ЗА ПОМОЩЬ!!!

Зосима · 25.05.2013, 15:18

Ну не знаю, почему отстой, мне нравится

Кстать, вместо real(I2) можно попробовать брать модуль abs(I2)

@BitMarshal · 27.05.2013, 01:06 **[ТС]**

Появилась ещё идея. Вот есть у меня трёхмерная матрица I(i,j,q). Как мне выделить из неё слой z = d (например z = 15), т.е. получить двухмерную матрицу Iz(i,j), полученную из слоя z = d, чтобы построить график ( X,Y,real(Iz) )? Мне это нужно, чтобы 30 раз не считать одну и ту же матрицу, а сразу строить графики при разных z.
P.S.: а можно ли как-то выделить в программе этот параметр d, чтобы иметь возможность сразу строить 2 графика при разных z (т.е. иметь параметры d1, d2)?

Добавлено через 22 часа 50 минут
Понимаю, что можно обойтись и без этого, но так мне удастся сохранить уйму времени.

Зосима · 27.05.2013, 13:36

Конечно можно

кстать, ты смотрел функцию slice?
Пусть мы уже расcчитали матрицу I(i,j,q) (ее можно один раз посчитать, сохранить в mat-файл, а потом просто загружать!

)
Пусть z = 0:n, d1 = n/3, d2 = 22.8, тогда:

Matlab M
1
2
3
4
5
6
% находим индекс эелмента z, наиболее близкого к d
[musor, id1] = min( abs(z-d1) ); 
[musor, id2] = min( abs(z-d2) );
% двумерная матрица с нужным слоем
In1 = I(:,:,id1); 
In2 = I(:,:,id2);

@BitMarshal · 29.05.2013, 01:03 **[ТС]**

Уважаемый Зосима!

Спасибо за Вашу помощь в написании программы! Она работает, поэтому текущая работа направлена на исправление формулы.
Функцию slice я смотрел, но она, к сожалению, не подходит для выполнения моей задачи.
В связи с задачами дипломной работы у меня возникли такие вопросы:

1. Что такое "musor", представленный в вашей последней программе? Это специальное название встроенной функции или прикол? Можете поподробнее рассказать о строке

Matlab M
1
[musor, id1] = min( abs(z-d1) );

Выражение справа означает минимум модуля разности z и d1 (только что означает "минимум" в данном контексте? Округление вниз, верно?), но что означает матрица слева?

2 (более важно). Мне необходимо найти поэлементную сумму матриц. Сначала двумерных, потом трёхмерных.
Подробнее: пусть есть одинаковые матрицы I1(i,j) и I2(i,j) с разными значениями. Необходимо получить матрицу I3(i,j), чтобы каждый её элемент был суммой соответствующих элементов матриц I1 и I2. Т.е

Matlab M
1
2
3
4
I1 = [1,2,3;4,5,6;7,8,9];
I2 = [10,11,12;13,14,15;16,17,18];
% Поэлементное суммирование матриц
I3 = [11,13,15;17,19,21;23,25,27]

Потом то же самое проделать с трёхмерными матрицами.
Это делается так, верно?

Matlab M
1
I3 = I1 + I2

С двухмерными матрицами получилось, но команда кажется мне слишком простой, чтобы быть правильной. Хотя..
Кстати, не подскажете, как можно задать трёхмерную матрицу в командном окне? Например размером 3х3х3?

3. Как сделать сетку более редкой? Пытался воспользоваться способом, описанным в этой статье, но у меня не получилось. Matlab не хотел делать сетку через каждые 5 единиц (10 столбцов матрицы), к тому же рисовал линии только по x.
Вот моя программа:

Matlab M
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
clc
clear all
 
% Начальные условия
a = 2; R = 50; L = 50;
k = 1; R0 = 4;
x0 = 15; y0 = 15; z0 = 15;
n1 = -40; n2 = 90; % Размеры сетки в плоскости Oxy
m = 20; % Количество шагов интегрирования
l = 10; % Количество линий уровня
z = 15; % Выбор слоя по координате z
 
% Расчёт
x = n1:0.5:n2; y = n1:0.5:n2; % Пределы изменения координат x и y
I = zeros(length(x),length(y));
for i = 1:length(x)
  for j = 1:length(y)
    [xx yy zz] = meshgrid( linspace(0,R,m), linspace(0,R,m), linspace(0,R,m) );
    r = sqrt( (x(i)-xx).^2 + (y(j)-yy).^2 + (z-zz).^2 );
    fr = z0 + sqrt(R0^2 - (xx-x0).^2 - (yy-y0).^2);
    f = (1 - cos(asin(a./r))).* fr * k * R^2 / (4 * (L + R)^2);
    f( isnan(f)|isinf(f) ) = 0;
    Fdz = trapz(linspace(0,R,m), f, 3);
    Fdydz = trapz(yy(:,1,1), Fdz);
    I(i,j) = trapz(xx(1,:,1), Fdydz);
  end
end
[X,Y] = meshgrid(x,y);
 
% Графики
figure('Color','w');
 
% Трехмерный график
subplot (1,2,2); surf(X,Y,real(I));
grid on; colormap jet;
shading faceted; colorbar;
alpha 0.8;
xlabel ('x'); ylabel ('y'); zlabel ('I');
title ('Свёртка');
xlim([n1 n2]); ylim([n1 n2]);
view (-17,19);
 
% График линий уровня
subplot (1,2,1); contourf (X,Y,real(I),l,'k');
xlabel ('x'); ylabel ('y');
title ('Свёртка (линии уровня)');
axis square;

А вот рисунок:

Зосима · 29.05.2013, 08:32

Пока еду, отвечу что знаю

1 musor - это просто имя переменной, которая не используется. Если ты заглянешь в хелп функции min, то увидишь, что когда принимается два параметра, то первый - это сам минимум, а второй - номер этого минимального элемента, который нам и нужен!
Собственно задача стоит в нахождении номера строки, где наиболее близко z=d (ведь значение d может и не быть в массиве z!).
Т.е. другими словами мы ищем номер элемента, где будет минимум разности z-d, только не просто разности, а модуля! Поэтому мы и выходим на такую формулу

понимаешь?

2 несмотря на кажущуюся простоту, операция "+" делает поэлементное сложение матриц произвольной размерности, а функция sum считает сумму элементов произвольной матрицы вдоль произвольной размерности.
Это можно считать фишкой матлаба - МАТричной ЛАБоратории - заточенной под работу с матрицами!
Можешь соорудить тестовую программку с двойным или тройным циклом
I3(i,j) = I2(i,j)+I1(i,j);
И убедиться, что цикл считает также, только во сто крат дольше

Если память не изменяет, то трехмерную матрицу без цикла можно задать так:
zeros(3,3,3), rand(3,3,3)
Но надо проверять (может rand([3,3,3]) )

3 тут пока не могу толком глянуть, разобраться. Хотя рисунок мне нравится

единственное линии уровня в surf можно былобы убрать.
После команды surf пропиши:
shading interp

*можно на "ты"

@BitMarshal · 04.06.2013, 18:30 **[ТС]**

Снова здравствуйте, Зосима!

С редкой сеткой я разобрался, пределы интегрирования также нашли похожие на правильные.
Но осталась такая проблема. Дело в том, что ваша программа в результирующей матрице перепутала местами значения x и y, т.е., она даёт матрицу I(y,x), а надо I(x,y). как это исправить?

@BitMarshal · 04.06.2013, 18:34 **[ТС]**

Вот какой график она даёт:

А надо

Решить эту проблему мне удалось, используя в качестве результата вычислений интегралов матрицу I(j,i), а не I(i,j). Но почему?

@BitMarshal · 04.06.2013, 18:40 **[ТС]**

Вот текст моей программы:

Matlab M
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
clc
clear all
 
% Исходные данные
a = 0.15; R = 50; L = 3.8; k = 1;
R1 = 4; x1 = 10; y1 = 20; z1 = 15;
R2 = 5; x2 = 30; y2 = 20; z2 = 15;
n1 = 0; n2 = 40;
m = 20; l = 10; z = 15;
t = 1; % Шаг координатной сетки
 
% Расчёт
x = n1:t:n2; y = n1:t:n2;
I1 = zeros(length(x),length(y));
for i1 = 1:length(x)
  for j1 = 1:length(y)
    [xx1 yy1 zz1] = meshgrid( linspace(x1-R1,x1+R1,m), linspace(y1-R1,y1+R1,m), linspace(z1-R1,z1+R1,m) );
    r1 = sqrt( (x(i1)-xx1).^2 + (y(j1)-yy1).^2 + (z-zz1).^2 );
    fr1 = z1 + sqrt(R1^2 - (xx1-x1).^2 - (yy1-y1).^2);
    f1 = (1 - cos(asin(a./r1))).* (1 ./ (L + R - r1).^2 +...
      1 ./ (L + R + r1).^2).* fr1 * k * R^2 / 8;
    f1( isnan(f1)|isinf(f1) ) = 0;
    Fdz1 = trapz(linspace(z1-R1,z1+R1,m), f1, 3);
    Fdydz1 = trapz(yy1(:,1,1), Fdz1);
    I1(j1,i1) = trapz(xx1(1,:,1), Fdydz1);
  end
end
I2 = zeros(length(x),length(y));
for i2 = 1:length(x)
  for j2 = 1:length(y)
    [xx2 yy2 zz2] = meshgrid( linspace(x2-R2,x2+R2,m), linspace(y2-R2,y2+R2,m), linspace(z2-R2,z2+R2,m) );
    r2 = sqrt( (x(i2)-xx2).^2 + (y(j2)-yy2).^2 + (z-zz2).^2 );
    fr2 = z2 + sqrt(R2^2 - (xx2-x2).^2 - (yy2-y2).^2);
    f2 = (1 - cos(asin(a./r2))).* (1 ./ (L + R - r2).^2 +...
      1 ./ (L + R + r2).^2).* fr2 * k * R^2 / 8;
    f2( isnan(f2)|isinf(f2) ) = 0;
    Fdz2 = trapz(linspace(z2-R2,z2+R2,m), f2, 3);
    Fdydz2 = trapz(yy2(:,1,1), Fdz2);
    I2(j2,i2) = trapz(xx2(1,:,1), Fdydz2);
  end
end
I3 = I1 + I2;
[X,Y] = meshgrid(x,y);
 
% Графики
figure('Color','w');
% Трехмерный график
subplot (1,2,2); surf(X,Y,real(I3) ,'EdgeColor','none');
grid on; colormap jet; alpha 1;
hold on; [ms,ns] = size(X);
is = 1:1/t:ns; plot3(X(:,is),Y(:,is),real(I3(:,is)),'k');
js = 1:1/t:ms; plot3(X(js,:)',Y(js,:)',real(I3(js,:))','k');
hold off;
xlabel ('x'); ylabel ('y'); zlabel ('I');
title ('Свёртка');
xlim([n1 n2]); ylim([n1 n2]); view (-17,19);
axis square;
% График линий уровня
subplot (1,2,1); contourf (X,Y,real(I3),l,'k');
xlabel ('x'); ylabel ('y');
title ('Свёртка (линии уровня)'); axis square;

Добавлено через 3 минуты
И ещё. Если внимательно посмотреть на правую возвышенность левого графика, можно увидеть круг в самом центре (по линии). Но чем дальше от центра холма, тем больше он становится похожим на квадрат. Потом снова круг (на границе между голубым и синим). Почему так происходит? Мне надо, что все линии вокруг центра второго холма должны быть кругами. Левый холм никак не влияет га правый, так как они рассчитывались независимо. Может, проблема в математике, а может.. Я не знаю.

Зосима · 04.06.2013, 18:50

*Исправь строки 14 и 28, вначале должно идти length(y)

Если записать I(i,j), то i - номер строки, а j - номер столбца, ОДНАКО (внимание!) номер строки определяет положение точки по вертикали! А номер столбца - по горизонтали!

(если непонятно, возьми листок в клетку, нарисуй матрицу 5х5, пронумеруй строки, столбцы (левый верхний угол - 1,1), понаходи координаты разных клеточек, закрой, забудь и жди просветления!

)

А мы привыкли к декартовой системе, где по вертикали меняется Y, а по горизонтали меняется X

Этот момент я упустил.

@BitMarshal · 04.06.2013, 20:16 **[ТС]**

Неа, не работает!
Программа:

Matlab M
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
clc
clear all
 
% Исходные данные
a = 0.15; R = 50; L = 3.8; k = 1;
R1 = 4; x1 = 10; y1 = 20; z1 = 15;
R2 = 4.5; x2 = 30; y2 = 20; z2 = 15;
n1 = 0; n2 = 40;
m = 20; l = 10; z = 15;
t = 1; % Шаг координатной сетки
 
% Расчёт
x = n1:t:n2; y = n1:t:n2;
I1 = zeros(length(y),length(x));
for i1 = 1:length(x)
  for j1 = 1:length(y)
    [xx1 yy1 zz1] = meshgrid( linspace(x1-R1,x1+R1,m), linspace(y1-R1,y1+R1,m), linspace(z1-R1,z1+R1,m) );
    r1 = sqrt( (x(i1)-xx1).^2 + (y(j1)-yy1).^2 + (z-zz1).^2 );
    fr1 = z1 + sqrt(R1^2 - (xx1-x1).^2 - (yy1-y1).^2);
    f1 = (1 - cos(asin(a./r1))).* (1 ./ (L + R - r1).^2 +...
      1 ./ (L + R + r1).^2).* fr1 * k * R^2 / 8;
    f1( isnan(f1)|isinf(f1) ) = 0;
    Fdz1 = trapz(linspace(z1-R1,z1+R1,m), f1, 3);
    Fdydz1 = trapz(yy1(:,1,1), Fdz1);
    I1(i1,j1) = trapz(xx1(1,:,1), Fdydz1);
  end
end
I2 = zeros(length(y),length(x));
for i2 = 1:length(x)
  for j2 = 1:length(y)
    [xx2 yy2 zz2] = meshgrid( linspace(x2-R2,x2+R2,m), linspace(y2-R2,y2+R2,m), linspace(z2-R2,z2+R2,m) );
    r2 = sqrt( (x(i2)-xx2).^2 + (y(j2)-yy2).^2 + (z-zz2).^2 );
    fr2 = z2 + sqrt(R2^2 - (xx2-x2).^2 - (yy2-y2).^2);
    f2 = (1 - cos(asin(a./r2))).* (1 ./ (L + R - r2).^2 +...
      1 ./ (L + R + r2).^2).* fr2 * k * R^2 / 8;
    f2( isnan(f2)|isinf(f2) ) = 0;
    Fdz2 = trapz(linspace(z2-R2,z2+R2,m), f2, 3);
    Fdydz2 = trapz(yy2(:,1,1), Fdz2);
    I2(i2,j2) = trapz(xx2(1,:,1), Fdydz2);
  end
end
I3 = I1 + I2;
[X,Y] = meshgrid(x,y);
 
% Графики
figure('Color','w');
% Трехмерный график
subplot (1,2,2); surf(X,Y,real(I3) ,'EdgeColor','none');
grid on; colormap jet; alpha 1;
hold on; [ms,ns] = size(X);
is = 1:1/t:ns; plot3(X(:,is),Y(:,is),real(I3(:,is)),'k');
js = 1:1/t:ms; plot3(X(js,:)',Y(js,:)',real(I3(js,:))','k');
hold off;
xlabel ('x'); ylabel ('y'); zlabel ('I');
title ('Свёртка');
xlim([n1 n2]); ylim([n1 n2]); view (-17,19);
axis square;
% График линий уровня
subplot (1,2,1); contourf (X,Y,real(I3),l,'k');
xlabel ('x'); ylabel ('y');
title ('Свёртка (линии уровня)'); axis square;

График:

@BitMarshal · 04.06.2013, 20:27 **[ТС]**

Замена

Matlab M
1
2
for i1 = 1:length(x)
  for j1 = 1:length(y)

на

Matlab M
1
2
for i1 = 1:length(y)
  for j1 = 1:length(x)

также ничего не меняет.

Добавлено через 9 минут
Ещё вопрос: если вывести строку

Matlab M
1
[xx yy zz] = meshgrid( linspace(x0-R0,x0+R0,m), linspace(y0-R0,y0+R0,m), linspace(z0-R0,z0+R0,m) );

из-под цикла, это приведёт к ускорению вычислений?

Зосима · 04.06.2013, 20:32

А так - рисует

Matlab M
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
clear, clc
 
clc
clear all
 
% Исходные данные
a = 0.15; R = 50; L = 3.8; k = 1;
R1 = 4; x1 = 10; y1 = 20; z1 = 15;
R2 = 4.5; x2 = 30; y2 = 20; z2 = 15;
n1 = 0; n2 = 40;
m = 20; l = 10; z = 15;
t = 1; % Шаг координатной сетки
 
% Расчёт
x = n1:t:n2; y = n1:t:n2;
I1 = zeros(length(y),length(x));
for i1 = 1:length(y)
  for j1 = 1:length(x)
    [xx1 yy1 zz1] = meshgrid( linspace(x1-R1,x1+R1,m), linspace(y1-R1,y1+R1,m), linspace(z1-R1,z1+R1,m) );
    r1 = sqrt( (x(j1)-xx1).^2 + (y(i1)-yy1).^2 + (z-zz1).^2 );
    fr1 = z1 + sqrt(R1^2 - (xx1-x1).^2 - (yy1-y1).^2);
    f1 = (1 - cos(asin(a./r1))).* (1 ./ (L + R - r1).^2 +...
      1 ./ (L + R + r1).^2).* fr1 * k * R^2 / 8;
    f1( isnan(f1)|isinf(f1) ) = 0;
    Fdz1 = trapz(linspace(z1-R1,z1+R1,m), f1, 3);
    Fdydz1 = trapz(yy1(:,1,1), Fdz1);
    I1(i1,j1) = trapz(xx1(1,:,1), Fdydz1);
  end
end
I2 = zeros(length(y),length(x));
for i2 = 1:length(y)
  for j2 = 1:length(x)
    [xx2 yy2 zz2] = meshgrid( linspace(x2-R2,x2+R2,m), linspace(y2-R2,y2+R2,m), linspace(z2-R2,z2+R2,m) );
    r2 = sqrt( (x(j2)-xx2).^2 + (y(i2)-yy2).^2 + (z-zz2).^2 );
    fr2 = z2 + sqrt(R2^2 - (xx2-x2).^2 - (yy2-y2).^2);
    f2 = (1 - cos(asin(a./r2))).* (1 ./ (L + R - r2).^2 +...
      1 ./ (L + R + r2).^2).* fr2 * k * R^2 / 8;
    f2( isnan(f2)|isinf(f2) ) = 0;
    Fdz2 = trapz(linspace(z2-R2,z2+R2,m), f2, 3);
    Fdydz2 = trapz(yy2(:,1,1), Fdz2);
    I2(i2,j2) = trapz(xx2(1,:,1), Fdydz2);
  end
end
I3 = I1 + I2;
[X,Y] = meshgrid(x,y);
 
% Графики
figure('Color','w');
% Трехмерный график
subplot (1,2,2); surf(X,Y,real(I3) ,'EdgeColor','none');
grid on; colormap jet; alpha 1;
hold on; [ms,ns] = size(X);
is = 1:1/t:ns; plot3(X(:,is),Y(:,is),real(I3(:,is)),'k');
js = 1:1/t:ms; plot3(X(js,:)',Y(js,:)',real(I3(js,:))','k');
hold off;
xlabel ('x'); ylabel ('y'); zlabel ('I');
title ('Свёртка');
xlim([n1 n2]); ylim([n1 n2]); view (-17,19);
axis square;
% График линий уровня
subplot (1,2,1); contourf (X,Y,real(I3),l,'k');
xlabel ('x'); ylabel ('y');
title ('Свёртка (линии уровня)'); axis square;

*заменил x(i) y(j) на x(j) y(i) ( ну и это оставил for i1 = 1:length(y), for j1 = 1:length(x) )

@BitMarshal · 04.06.2013, 20:43 **[ТС]**

Вот смотрите. У вас i - это индекс точки по координате y, j - по координате x. Поэтому матрица I(j,i) = I(y,x), как и в моём решении, когда я просто использовал I(j,i)! Вам не кажется, что вы в результате сделали то же самое, но НАМНОГО дольше?

Зосима · 04.06.2013, 21:03

Сообщение от BitMarshal

Ещё вопрос: если вывести строку

если у тебя пределы не зависят от текущих x,y то можно!

Добавлено через 15 минут

Сообщение от BitMarshal

вы в результате сделали то же самое, но НАМНОГО дольше?

может быть.... я после домашней лесопилки уже мало что соображаю

@BitMarshal · 05.06.2013, 11:12 **[ТС]**

Всё! Мне кажется, я добился правильной программы. Как и обещал, выкладываю её полный текст, чтобы все знали, как правильно программировать трёхмерную свёртку двух функций (имея в результате 4-мерную функцию, зависящую от x, y, z), вычисляемую методом трапеций, с выводом результата в виде графика значений, взятых в определённом слое z.

Matlab M
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
clc
clear all
 
% Исходные данные
    % На этой строке можно добавить дополнительные параметры подынтегральной функции
n1 = -20; n2 = 20; % Границы координатной сетки
m = 50; % Количество шагов интегрирования
t = 0.1; % Шаг координатной сетки
z = 1; % Выбор слоя трёхмерного пространства
 
% Расчёт
x = n1:t:n2; y = n1:t:n2;
I = zeros(length(y),length(x));
for i = 1:length(x)
  for j = 1:length(y)
    [xx yy zz] = meshgrid( linspace(A1,A2,m), linspace(B1,B2,m), linspace(C1,C2,m) );
    f = FFF;
    f( isnan(f)|isinf(f) ) = 0;
    Fdz = trapz(linspace(C1,C2,m), f, 3);
    Fdydz = trapz(yy(:,1,1), Fdz);
    F(j,i) = trapz(xx(1,:,1), Fdydz);
  end
end
[X,Y] = meshgrid(x,y);
 
% Трехмерный график
surf(X,Y,real(F));
grid on; shading interp;
xlabel ('x'); ylabel ('y'); zlabel ('F');
title ('Название графика');
xlim([n1 n2]); ylim([n1 n2]);
axis square;

Примечание:
A1, A2, B1, B2, C1, C2 - пределы интегрирования,
FFF - подынтегральная функция (в виде выражения, зависящего от x, y, z + дополнительных параметров),
F - результат вычисления функции свёртки.

Уважаемый Зосима, спасибо Вам большое за всё!
Надеюсь, проблем у меня больше не возникнет.

Добавлено через 10 часов 42 минуты
Уважаемый Зосима! Я вынес матрицу [xx yy zz] из-под цикла в своей дипломной работе, но последняя программа, выложенная мной в предыдущем сообщении, представляла собой общий вид программы такого рода. Как Вы справедливо заметили, упомянутую матрицу можно вынести, если только пределы интегрирования не зависят от переменных x, y, z. Но в общем случае это необходимо учитывать, поэтому я не стал выносить матрицу.

@BitMarshal · 08.06.2013, 22:38 **[ТС]**

Доброго времени суток! Снова..

Теперь у меня вопросы чисто по оформлению.

1. Как заставить Matlab отображать в качестве подписей осей графика Юникод-символы? Т.е. вместо x и y мне надо написать φ и θ соответственно.

2. Как нарисовать цветовую шкалу сразу под двумя графиками, во всю длину окна? Команда

Matlab M
1
colorbar('SouthOutside')

приводит к графику

но мне нужно, чтобы он был посредине окна сразу под двумя графиками. Я, конечно. понимаю, что это можно сделать вручную в окне Figure, но мне нужно прописать это программно.

Новые блоги и статьи Все статьи Все блоги /
Новый ноутбук volvo 07.12.2025 Всем привет. По скидке в "черную пятницу" взял себе новый ноутбук Lenovo ThinkBook 16 G7 на Амазоне: Ryzen 5 7533HS 64 Gb DDR5 1Tb NVMe 16" Full HD Display Win11 Pro	Музыка, написанная Искусственным Интеллектом volvo 04.12.2025 Всем привет. Некоторое время назад меня заинтересовало, что уже умеет ИИ в плане написания музыки для песен, и, собственно, исполнения этих самых песен. Стихов у нас много, уже вышли 4 книги, еще 3. . .	От async/await к виртуальным потокам в Python IndentationError 23.11.2025 Армин Ронахер поставил под сомнение async/ await. Создатель Flask заявляет: цветные функции - провал, виртуальные потоки - решение. Не threading-динозавры, а новое поколение лёгких потоков. Откат?. . .	Поиск "дружественных имён" СОМ портов Argus19 22.11.2025 Поиск "дружественных имён" СОМ портов На странице: https:/ / norseev. ru/ 2018/ 01/ 04/ comportlist_windows/ нашёл схожую тему. Там приведён код на С++, который показывает только имена СОМ портов, типа,. . .	Сколько Государство потратило денег на меня, обеспечивая инсулином. Programma_Boinc 20.11.2025 Сколько Государство потратило денег на меня, обеспечивая инсулином. Вот решила сделать интересный приблизительный подсчет, сколько государство потратило на меня денег на покупку инсулинов. . . .
Ломающие изменения в C#.NStar Alpha Etyuhibosecyu 20.11.2025 Уже можно не только тестировать, но и пользоваться C#. NStar - писать оконные приложения, содержащие надписи, кнопки, текстовые поля и даже изображения, например, моя игра "Три в ряд" написана на этом. . .	Мысли в слух kumehtar 18.11.2025 Кстати, совсем недавно имел разговор на тему медитаций с людьми. И обнаружил, что они вообще не понимают что такое медитация и зачем она нужна. Самые базовые вещи. Для них это - когда просто люди. . .	Создание Single Page Application на фреймах krapotkin 16.11.2025 Статья исключительно для начинающих. Подходы оригинальностью не блещут. В век Веб все очень привыкли к дизайну Single-Page-Application . Быстренько разберем подход "на фреймах". Мы делаем одну. . .	Фото: Daniel Greenwood kumehtar 13.11.2025	Расскажи мне о Мире, бродяга kumehtar 12.11.2025 — Расскажи мне о Мире, бродяга, Ты же видел моря и метели. Как сменялись короны и стяги, Как эпохи стрелою летели. - Этот мир — это крылья и горы, Снег и пламя, любовь и тревоги, И бескрайние. . .

@BitMarshal 18 / 18 / 0 Регистрация: 15.05.2013 Сообщений: 73

	Построение графика трёхмерной свёртки двух функций в Matlab 15.05.2013, 22:43. Показов 10803. Ответов 71 Метки нет (Все метки) Доброго времени суток! Пишу первый раз, поэтому прошу не судить меня строго за плохое оформление. Задача связана с дипломной работой, поэтому решение необходимо срочно. А проблема звучит так: необходимо построить на одном рисунке два графика (линий уровня и трёхмерный) трёхмерной свёртки двух функций: функции источника единичной интенсивности f(x,y,z) и функции плотности распределения вещества g(x,y,z). Функция источника выглядит так: f(x,y,z) = [1 – cos (arcsin (a/r))] * [1 / (L + R – r)² + 1 / (L + R + r)²] * k * R² / 8, где r = √(x² + y² + z²). Функция плотности (имеет вид сферы): g(x,y,z) = z₀ + √(R₀² – (x – x₀)² – (y – y₀)²), где x₀, y₀, z₀ — настраиваемые координаты центра области распределения. Свёртка этих функций должна иметь такой вид: I (x,y,z) = ∫∫∫ f(x–x', y–y', z–z') * g(x', y', z') dx' dy' dz'. В случае моей задачи: I (x,y,z) = ∫∫∫ [1 – cos(arcsin(a / √((x–x')² + (y–y')² + (z–z')²)))] * [1 / (L + R – √((x–x')² + (y–y')² + (z–z')²))² + 1 / (L + R + √((x–x')² + (y–y')² + (z–z')²))²] * √(R₀² – (x'–x₀)² – (y'–y₀)²) * k * R² / 8 dx' dy' dz'. Пределы интегрирования — (0, x), (0, y) и (0, z) соответственно. Проведя численное интегрирование, необходимо построить график поверхности от (x,y,I). Координаты x,y,z должны менятся от -25 до 25 каждая. В начале программы также должен быть блок вводимых параметров: a = 0.01 (с клавиатуры), R = 50, L = 50, R₀ = 4, k = 1. Затем необходимо сделать всё то же самое, но уже в качестве функции плотности будет выступать свёртка двух сфер, т.е. ∫∫ [z₁ + √(R₁² – (x'–x₁)² – (y'–y₁)²)] * [z₂ + √(R₁² – (x–x₂–x')² – (y–y₂–y')²)] dx' dy'. Проблема заключается в том, что Matlab отказывается проводить интегрирование с помощью команды triplequad (явно интеграл посчитать не удаётся) с переменными пределами. Прошу мне помочь, так как для меня это очень важно. Matlabэом пользуюсь всего месяц и много чего не знаю. Версия моей программы — R2012b (8.00). 0

Зосима 5245 / 3573 / 379 Регистрация: 02.04.2012 Сообщений: 6,477 Записей в блоге: 18
	24.05.2013, 17:09
	так, хорошо, так что осталось делать? чтот я нить потерял... 0

Зосима 5245 / 3573 / 379 Регистрация: 02.04.2012 Сообщений: 6,477 Записей в блоге: 18
	25.05.2013, 13:41
	Ага кажется дошло ты правильно рассуждаешь. Попробуй: I2 = sum(I, 3); *здесь 3 - номер размерности, по которой идет суммирование. Поидее, если я не напортачил, то I2 должно получиться двумерной матричкой. Проверяй 1

@BitMarshal 18 / 18 / 0 Регистрация: 15.05.2013 Сообщений: 73
	25.05.2013, 14:41 [ТС]
	Да, ваш способ работает, I(i,j,q) стала I2(i,j). Большое спасибо! График, правда, всё равно отстой, но теперь я точно уверен, что ошибка в формуле. БЛАГОДАРЮ ЗА ПОМОЩЬ!!! 0

Зосима 5245 / 3573 / 379 Регистрация: 02.04.2012 Сообщений: 6,477 Записей в блоге: 18
	25.05.2013, 15:18
	Ну не знаю, почему отстой, мне нравится Кстать, вместо real(I2) можно попробовать брать модуль abs(I2) 0

@BitMarshal 18 / 18 / 0 Регистрация: 15.05.2013 Сообщений: 73
	27.05.2013, 01:06 [ТС]
	Появилась ещё идея. Вот есть у меня трёхмерная матрица I(i,j,q). Как мне выделить из неё слой z = d (например z = 15), т.е. получить двухмерную матрицу Iz(i,j), полученную из слоя z = d, чтобы построить график ( X,Y,real(Iz) )? Мне это нужно, чтобы 30 раз не считать одну и ту же матрицу, а сразу строить графики при разных z. P.S.: а можно ли как-то выделить в программе этот параметр d, чтобы иметь возможность сразу строить 2 графика при разных z (т.е. иметь параметры d1, d2)? Добавлено через 22 часа 50 минут Понимаю, что можно обойтись и без этого, но так мне удастся сохранить уйму времени. 0

Зосима 5245 / 3573 / 379 Регистрация: 02.04.2012 Сообщений: 6,477 Записей в блоге: 18
	29.05.2013, 08:32
	Пока еду, отвечу что знаю 1 musor - это просто имя переменной, которая не используется. Если ты заглянешь в хелп функции min, то увидишь, что когда принимается два параметра, то первый - это сам минимум, а второй - номер этого минимального элемента, который нам и нужен! Собственно задача стоит в нахождении номера строки, где наиболее близко z=d (ведь значение d может и не быть в массиве z!). Т.е. другими словами мы ищем номер элемента, где будет минимум разности z-d, только не просто разности, а модуля! Поэтому мы и выходим на такую формулу понимаешь? 2 несмотря на кажущуюся простоту, операция "+" делает поэлементное сложение матриц произвольной размерности, а функция sum считает сумму элементов произвольной матрицы вдоль произвольной размерности. Это можно считать фишкой матлаба - МАТричной ЛАБоратории - заточенной под работу с матрицами! Можешь соорудить тестовую программку с двойным или тройным циклом I3(i,j) = I2(i,j)+I1(i,j); И убедиться, что цикл считает также, только во сто крат дольше Если память не изменяет, то трехмерную матрицу без цикла можно задать так: zeros(3,3,3), rand(3,3,3) Но надо проверять (может rand([3,3,3]) ) 3 тут пока не могу толком глянуть, разобраться. Хотя рисунок мне нравится единственное линии уровня в surf можно былобы убрать. После команды surf пропиши: shading interp *можно на "ты" 0

@BitMarshal 18 / 18 / 0 Регистрация: 15.05.2013 Сообщений: 73
	04.06.2013, 18:30 [ТС]
	Снова здравствуйте, Зосима! С редкой сеткой я разобрался, пределы интегрирования также нашли похожие на правильные. Но осталась такая проблема. Дело в том, что ваша программа в результирующей матрице перепутала местами значения x и y, т.е., она даёт матрицу I(y,x), а надо I(x,y). как это исправить? 0

@BitMarshal 18 / 18 / 0 Регистрация: 15.05.2013 Сообщений: 73
	04.06.2013, 18:34 [ТС]
	Вот какой график она даёт: А надо Решить эту проблему мне удалось, используя в качестве результата вычислений интегралов матрицу I(j,i), а не I(i,j). Но почему? 0

Зосима 5245 / 3573 / 379 Регистрация: 02.04.2012 Сообщений: 6,477 Записей в блоге: 18
	04.06.2013, 18:50
	Исправь строки 14 и 28, вначале должно идти length(y) Если записать I(i,j), то i - номер строки, а j - номер столбца, ОДНАКО (внимание!) номер строки определяет положение точки по вертикали! А номер столбца - по горизонтали*! (если непонятно, возьми листок в клетку, нарисуй матрицу 5х5, пронумеруй строки, столбцы (левый верхний угол - 1,1), понаходи координаты разных клеточек, закрой, забудь и жди просветления! ) А мы привыкли к декартовой системе, где по вертикали меняется Y, а по горизонтали меняется X Этот момент я упустил. 0

@BitMarshal 18 / 18 / 0 Регистрация: 15.05.2013 Сообщений: 73
	04.06.2013, 20:43 [ТС]
	Вот смотрите. У вас i - это индекс точки по координате y, j - по координате x. Поэтому матрица I(j,i) = I(y,x), как и в моём решении, когда я просто использовал I(j,i)! Вам не кажется, что вы в результате сделали то же самое, но НАМНОГО дольше? 0