Нужна помощь в построении этой функции в Matlab

Недавно я начал использовать Matlab и пытаюсь построить мнимую часть функции. Я вижу, что где-то делаю ошибку, потому что у меня есть график, показывающий, что мне нужно получить, а я получаю что-то другое. Вот ссылка на график, который мне нужно получить, который я получаю, и функцию, которую я рисую:

Искомые результаты, построенная функция и полученный график.

Я знаю, что эта функция имеет сингулярность на частоте около 270 Гц, и я не знаю, сможет ли «квадроцикл» решить интеграл в таком случае. Вы, ребята, наверное, знаете, но тета-функция внутри интеграла — это хевисайд, и я не знаю, может быть, это вызывает проблему? Что-то не так, что я здесь делаю? Вот код Matlab, который я написал:

clear all
clc

ff=1:10:1000;

K1=(2*3)*log(2*cosh(135/6))/pi;
theta1=ones(size(ff));
theta1(ff<270)=0;

I1=zeros(size(ff));

for n = 1:numel(ff)
    f = ff(n);
    I1(n)= K1/f + (f/pi)*quadgk(@(x)(sinh(x/3)/(cosh(135/3)+cosh(x/3))-theta1(n))./((f^2)-4*(x.^2)), 0, inf);
end

plot(ff,I1, 'b');
hold on
axis([0 1000 -0.3 1])

function matlab plot integral singular
person user2768562    schedule 11.09.2013    source источник
comment
Я посмотрел на формулу и вашу реализацию в Matlab, в этом нет ничего плохого (что, кстати, также предполагает сюжет). Так что я думаю, вам следует попробовать найти альтернативу quadgk, так как это, похоже, проблема. Я не могу помочь вам с этим, хотя.   -  person Fraukje    schedule 11.09.2013
comment
Я нашел это в документации Matlab quad (см. mathworks.nl/help/matlab /ref/quad.html): если функция сингулярна в точках внутри (a,b), запишите интеграл как сумму интегралов по подинтервалам с особыми точками в качестве конечных точек, вычислите их с помощью quadgk и добавьте результаты, достижения. Тем не менее, все еще кажется, что в вашей функции есть несколько особенностей...?   -  person Fraukje    schedule 11.09.2013
comment
ну, так как theta прерывистый, мы найдем этот шаг и в функции I1.   -  person fpe    schedule 11.09.2013

Ответы (2)


arrow_upward
2
arrow_downward

  • Во-первых, я изменил выражение для вашей функции Хевисайда, придав ему правильную форму.
  • Во-вторых, я явно ввел определения mu и T.
  • В-третьих, я разбил интеграл на 4 компонента следующим образом и оценил их по отдельности (в соответствии с тем, что предложил Фраукье).
  • В-четвертых, я использую quadl, хотя это проблема. второстепенного значения.
  • (изменить) Изменен диапазон ff

Вот код с изменениями:

fstep=1;
ff=[1:fstep:265 275:fstep:1000];

T = 3;
mu = 135;

df = 0.01;
xmax = 10;

K1=(2*T/pi)*log(2*cosh(mu/(2*T)));
theta1=ones(size(ff));
theta1(ff-2*mu<0)=0;

I1=zeros(size(ff));    
for n = 1:numel(ff)
    f = ff(n);
    sigm1 = @(x)  sinh(x/T)./((f^2-4*x.^2).*(cosh(mu/T)+cosh(x/T)));
    sigm2 = @(x)   -theta1(n)./(f^2-4*x.^2);
    I1(n) = K1/f  + (f/pi)*quadl(sigm1,0,f/2-df);      % term #1
%    I1(n) = I1(n) + (f/pi)*quadl(sigm1,f/2+df,xmax);   % term #2
%    I1(n) = I1(n) + (f/pi)*quadl(sigm2,0,f/2-df);     % term #3
%    I1(n) = I1(n) + (f/pi)*quadl(sigm2,f/2+df,xmax);   % term #4
end

Я решил разделить интегралы вокруг x=f/2, так как там явно есть сингулярность (деление на 0). Но дополнительные проблемы возникают для членов № 2 и № 4, то есть когда интегралы оцениваются для x>f/2, предположительно из-за всех тригонометрических членов.

Если вы сохраните только термины 1 и 3, вы получите что-то достаточно похожее на показанный вами сюжет:

введите здесь описание изображения

Однако вам, вероятно, следует более внимательно изучить свою функцию и посмотреть, что можно сделать для вычисления интегралов для x>f/2.

ИЗМЕНИТЬ

Я еще раз просмотрел код и переопределил вспомогательные интегралы:

I1=zeros(size(ff));
I2=zeros(size(ff));
I3=zeros(size(ff));

for n = 1:numel(ff)
    f = ff(n);
    sigm3 = @(x)  sinh(x/T)./((f^2-4*x.^2).*(cosh(mu/T)+cosh(x/T))) -theta1(n)./(f^2-4*x.^2);
    I1(n) = K1/f  + (f/pi)*quadl(sigm3,0,f/2-df); 
    I2(n) = (f/pi)*quadl(sigm3,f/2+df,10);
end
I3=I2;
I3(isnan(I3)) = 0;
I3 = I3 + I1;

Вот как теперь выглядит вывод:

введите здесь описание изображения

Зеленая линия представляет собой интеграл функции по 0<x<f/2 и выглядит нормально. Красная линия представляет собой интеграл по Inf>x>f/2 и явно не соответствует f=270. Синяя кривая представляет собой сумму (общий интеграл) без учета вклада NaN при интегрировании по Inf>x>f/2.

Мой вывод состоит в том, что может быть что-то не так с кривой, как вы ожидаете.

person Buck Thorn    schedule 11.09.2013

arrow_upward
1
arrow_downward

Пока я бы поступил так:

clc,clear

T = 3;
mu = 135;

f = 1E-04:.1:1000;

theta = ones(size(f));
theta(f < 270)= 0;

integrative = zeros(size(f));
for ii = 1:numel(f)
    ff = @(x) int_y(x, f(ii), theta(ii));
    integrative(ii) = quad(ff,0,2000);
end

Imm = ((2*T)./(pi*f)).*log(2*cosh(mu/(2*T))) + (f/pi).*integrative;

Imm1 = exp(interp1(log(f(1:2399)),log(Imm(1:2399)),log(f(2400):.001:f(2700)),'linear','extrap'));
Imm2 = exp(interp1(log(f(2985:end)),log(Imm(2985:end)),log(f(2701):.001:f(2984)),'linear','extrap'));


plot([(f(2400):.001:f(2700)) (f(2701):.001:f(2984))],[Imm1 Imm2])
hold on
axis([0 1000 -1.0 1])
plot(f,Imm,'g')
grid on
hold off

с

function rrr = int_y(x,f,theta)
T = 3;
mu = 135;
rrr = ( (sinh(x./T)./(cosh(mu/T) + cosh(x/T))) - theta ) ./ (f.^2 - 4.*(x.^2));
end

Я придумал такой сюжет:

введите здесь описание изображения

person fpe    schedule 11.09.2013