5

我最近开始使用 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])
4

2 回答 2

2
  • 首先，我将您的重载函数的表达式更改为我认为正确的形式。
  • 其次，我明确地介绍了 mu 和 T 的定义。
  • 第三，我将积分分解为 4 个组件，如下所示，并分别评估它们（按照 Fraukje 提出的思路）
  • 第四，我使用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。

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

于 2013-09-11T14:11:17.637 回答
1

到目前为止，我会这样进行：

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

我想出了这个情节：

在此处输入图像描述

于 2013-09-11T14:03:58.300 回答