Двойной интеграл в MATLAB

Question

Я пытаюсь найти мощность в апертуре от гауссова луча, где апертура смещена от центра луча. Решением является следующее уравнение ( ссылка ) (извините, здесь нет латекса):

Wz является константой вместе с a и r . Я не уверен, как я могу сделать что-то подобное с MATLAB. У кого-нибудь есть предложение? Я знаю, что есть функция dblquad(), но она предполагает, что пределы интеграции фиксированы и не зависят друг от друга.

Answer 1

Используя немного математической работы, вы можете уменьшить двойной интеграл до единственного (хотя и содержащего функцию ошибки), который должен быть легче вычислен численно в MATLAB:

(с оговоркой на ошибки; если возможно, проверьте расчеты самостоятельно.)

Answer 2

Оказывается, что более поздние версии MATLAB теперь имеют функцию quad2d(), которая выполняет двумерный интеграл по поверхности. Пример 2 на странице ссылок подробно описывает пример такого типа интеграции.

Мой код в итоге выглядел примерно так:

powerIntegral = @(x,y) 2/(pi*W^2)*exp(-2*((x - offsetDist).^2 + y.^2)/(W^2)); 
ymin = @(x) -sqrt(radius.^2 - x.^2);
ymax = @(x) sqrt(radius.^2 - x.^2);
powerRatioGaussian = quad2d(powerIntegral,-radius,radius,ymin,ymax);

Довольно изящно.Спасибо за помощь.

Answer 3

Я не уверен, но я думаю, что символическая панель инструментов может помочь вам здесь. Он подходит для такого рода проблем. Вы можете определить свои переменные как символьные переменные, используя команду syms , и вычислить интеграл символически. Затем вы можете назначить значения переменных и найти фактическое значение.

Отказ от ответственности : На самом деле я никогда не использовал его сам.

Answer 4

Вообще говоря, для числовой интеграции вы можете преобразовать интеграл с зависимыми граничными условиями в интеграл с независимыми границами, умножив его на функцию, которая равна 1, если вы находитесь внутри исходной границы, и 0, если вы находитесь снаружи. Тогда возьмите свои пределы, чтобы быть квадратом, который содержит ваши исходные условия. Другими словами, здесь вы умножите на

г (х, у) = ((х ^ 2 + у ^ 2) <а ^ 2) </p>

и ваши пределы будут -a

Вы должны быть немного осторожнее с предположениями о непрерывности в вашем методе интеграции, но вы должны быть в порядке, если что-то не очень странное. Вы всегда можете проверить это, изменив размер ячейки и убедившись, что вычисленное интегральное значение не меняется.

В этом конкретном случае вы также можете выполнить преобразование из декартовых в полярные координаты;

x = rcos (т)

y = rsin (t)

dxdy = rdrdt

Тогда ваши пределы интегрирования будут r от 0 до a и t от 0 до 2 * pi