需要帮助在 Matlab 中绘制此函数答案

【问题标题】：Need help plotting this function in Matlab需要帮助在 Matlab 中绘制此函数
【发布时间】：2013-09-15 10:31:00
【问题描述】：

我最近开始使用 Matlab，我正在尝试绘制函数的虚部。我可以看到我在某个地方犯了一个错误，因为我有图表显示我需要得到什么并且我得到了其他东西。这是我需要获取的图表的链接，我正在获取以及我正在绘制的函数：

我知道这个函数在 270 Hz 左右的频率上有一个奇点，我不知道“quadgk”是否可以解决积分问题。你们可能知道，但积分内的 theta 函数是一个重磅炸弹，我不知道这是否会导致问题？我在这里做错了什么吗？这是我写的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])

【问题讨论】：

我盯着公式和你在 Matlab 中的实现，似乎没有什么问题（顺便说一下，情节也暗示了这一点）。所以我认为你应该尝试为quadgk 寻找替代方案，因为这似乎是问题所在。不过我帮不了你。
我在 quad 的 Matlab 文档中发现了这一点（参见 mathworks.nl/help/matlab/ref/quad.html）：“如果函数在 (a,b) 内部的点处是奇异的，则将积分写为积分的总和以奇异点为端点的子区间，用quadgk计算它们，并将结果相加。”但是，您的函数中似乎仍然存在多个奇点......？
好吧，由于theta 是不连续的，我们会在I1 函数中找到该步骤。

标签： matlab function plot integral singular

【解决方案1】：

首先，我将您的 heaviside 函数的表达式更改为我认为正确的形式。
其次，我明确介绍了mu和T的定义。
第三，我将积分分解为 4 个组件，如下所示，并分别评估它们（按照 Fraukje 提出的思路）
第四，我使用quadl，虽然这很可能。次要的。
(edit) 更改了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 附近失败。蓝色曲线是在对Inf>x>f/2 进行积分时不包括NaN 贡献的总和（总积分）。

我的结论是，曲线可能有问题，正如您所期望的那样。

【讨论】：

【解决方案2】：

到目前为止，我会这样做：

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

我想出了这个情节：

【讨论】：