5

私は最近Matlabの使用を開始し、関数の虚部をプロットしようとしています. 取得する必要があるものを示すグラフがあり、別のものを取得しているため、どこかで間違いを犯していることがわかります。これは、私が取得する必要があるグラフと、プロットしている関数のリンクです。

欲しい結果、プロットされた関数、および得られたグラフ。

この関数には 270 Hz 付近の周波数で特異点があることはわかっていますが、「quadgk」で積分を解くことができるかどうかはわかりません。皆さんはおそらくご存知でしょうが、積分内のシータ関数はヘヴィサイドであり、それが問題を引き起こしているかどうかはわかりませんか? ここで何か間違ったことはありますか？これが私が書いた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 に答える