2

このコードから開始：

    clc, clear all, close all
tic

k1 = 0.01:0.1:100;
k2 = 0.01:0.1:100;
k3 = 0.01:0.1:100;

k = sqrt(k1.^2 + k2.^2 + k3.^2);

c = 1.476;
gamma = 3.9;

colors = {'cyan'};
Ek = (1.453*k.^4)./((1 + k.^2).^(17/6));
E = @(k) (1.453*k.^4)./((1 + k.^2).^(17/6));
E_int = zeros(1,numel(k1));
E_int(1) = 1.5;

for i = 2:numel(k)
    E_int(i) = E_int(i-1) - integral(E,k(i-1),k(i));
end

beta = c*gamma./(k.*sqrt(E_int));


F_11 = zeros(1,numel(k1));
F_22 = zeros(1,numel(k1));
F_33 = zeros(1,numel(k1));

count = 0;
for i = 1:numel(k1)
    count = count + 1;
    phi_11 = @(k2,k3) phi_11_new(k1,k2,k3,beta,i);
    phi_22 = @(k2,k3) phi_22_new(k1,k2,k3,beta,i);
    phi_33 = @(k2,k3) phi_33_new(k1,k2,k3,beta,i);
    F_11(count) = integral2(phi_11,-100,100,-100,100);
    F_22(count) = integral2(phi_22,-100,100,-100,100);
    F_33(count) = integral2(phi_33,-100,100,-100,100);
end

figure
hold on
plot(k1,F_11,'b')
plot(k1,F_22,'cyan')
plot(k1,F_33,'magenta')
hold off

どこ

function phi_11 = phi_11_new(k1,k2,k3,beta,i)
k = sqrt(k1(i).^2 + k2.^2 + k3.^2);
k30 = k3 + beta(i).*k1(i);
k0 = sqrt(k1(i).^2 + k2.^2 + k30.^2);
E_k0 = 1.453.*k0.^4./((1 + k0.^2).^(17/6));
C1 = (beta(i).*k1(i).^2).*(k1(i).^2 + k2.^2 - k3.*k30)./(k.^2.*(k1(i).^2 + k2.^2));
C2 = k2.*k0.^2./((k1(i).^2 + k2.^2).^(3/2)).*atan2((beta(i).*k1(i).*sqrt(k1(i).^2 + k2.^2)),(k0.^2 - k30.*k1(i).*beta(i)));
xhsi1 = C1 - k2./k1(i).*C2;
xhsi1_q = xhsi1.^2;
phi_11 = E_k0./(4.*pi.*k0.^4).*(k0.^2 - k1(i).^2 - 2.*k1(i).*k30.*xhsi1 + (k1(i).^2 + k2.^2).*xhsi1_q);
end

function phi_22 = phi_22_new(k1,k2,k3,beta,i)
k = sqrt(k1(i).^2 + k2.^2 + k3.^2);
k30 = k3 + beta(i).*k1(i);
k0 = sqrt(k1(i).^2 + k2.^2 + k30.^2);
E_k0 = 1.453.*k0.^4./((1 + k0.^2).^(17/6));
C1 = (beta(i).*k1(i).^2).*(k1(i).^2 + k2.^2 - k3.*k30)./(k.^2.*(k1(i).^2 + k2.^2));
C2 = k2.*k0.^2./((k1(i).^2 + k2.^2).^(3/2)).*atan2((beta(i).*k1(i).*sqrt(k1(i).^2 + k2.^2)),(k0.^2 - k30.*k1(i).*beta(i)));
xhsi2 = k2./k1(i).*C1 + C2;
xhsi2_q = xhsi2.^2;
phi_22 = E_k0./(4.*pi.*k0.^4).*(k0.^2 - k2.^2 - 2.*k2.*k30.*xhsi2 + (k1(i).^2 + k2.^2).*xhsi2_q);
end

function phi_33 = phi_33_new(k1,k2,k3,beta,i)
k = sqrt(k1(i).^2+k2.^2+k3.^2);
k30 = k3 + beta(i).*k1(i);
k0 = sqrt(k1(i).^2+k2.^2+k30.^2);
E_k0 = (1.453.*k0.^4./((1+k0.^2).^(17/6)));
phi_33 = (E_k0./(4*pi.*(k.^4))).*(k1(i).^2+k2.^2);
end

この手順は、研究から得られた他のいくつかと一致しない結果に私を導きます。一致させる必要のある結果を以下に示します。

ここに画像の説明を入力してください

私のはこれらのように見えますが

ここに画像の説明を入力してください

compwだけが理論結果とどのように一致するかを評価するのは非常に簡単です。したがって、この欠陥は、関数phi_11_new（およびphi_22_new）の外側のベータの定義にある可能性があると思います。

私が現在行っている方法ではなく、phi_11_new（およびphi_22_new）内でベータを計算する方法を提案する人はいますか？

よろしくお願いします。

よろしく、fpe

4

3 に答える 3

2

補間を改善して、小さな値で分解されないようにしました。また、値の対数を補間するようになったため、より正確な値が返されます。

これが現在のコードです。

function test15()

[k1,k2,k3] = deal(0.01:0.1:400);

k = sqrt(k1.^2 + k2.^2 + k3.^2);

c = 1.476;
gamma = 3.9;

Ek = (1.453*k.^4)./((1 + k.^2).^(17/6));
E_int = 1.5-cumtrapz(k,Ek);
beta = c*gamma./(k.*sqrt(E_int));

[F_11,F_22,F_33] = deal(zeros(1,numel(k1)));

k_vec = k;
beta_vec = beta;

kLim = 100;

for ii = 1:numel(k1)
    phi_11 = @(k2,k3) phi_11_new(k1(ii),k2,k3,k_vec,beta_vec);
    phi_22 = @(k2,k3) phi_22_new(k1(ii),k2,k3,k_vec,beta_vec);
    phi_33 = @(k2,k3) phi_33_new(k1(ii),k2,k3,k_vec,beta_vec);
    F_11(ii) = quad2d(phi_11,-kLim,kLim,-kLim,kLim);
    F_22(ii) = quad2d(phi_22,-kLim,kLim,-kLim,kLim);
    F_33(ii) = quad2d(phi_33,-kLim,kLim,-kLim,kLim);
end

figure
loglog(k1,F_11,'b')
hold on
loglog(k1,F_22,'cyan')
loglog(k1,F_33,'magenta')
hold off
grid on

end

function phi_11 = phi_11_new(k1,k2,k3,k_vec,beta_vec)
k = sqrt(k1^2 + k2.^2 + k3.^2);

log_beta_vec = interp1(log(k_vec),log(beta_vec),log(k(:)),'linear','extrap');
log_beta = reshape(log_beta_vec,size(k));
beta = exp(log_beta);

k30 = k3 + beta*k1;
k0 = sqrt(k1^2 + k2.^2 + k30.^2);
E_k0 = 1.453*k0.^4./((1 + k0.^2).^(17/6));
C1 = (beta*k1^2).*(k1^2 + k2.^2 - k3.*k30)./(k.^2.*(k1^2 + k2.^2));
C2 = k2.*k0.^2./((k1^2 + k2.^2).^(3/2)).*atan2((beta*k1.*sqrt(k1^2 + k2.^2)),(k0.^2 - k30*k1.*beta));
xhsi1 = C1 - (k2/k1).*C2;
xhsi1_q = xhsi1.^2;
phi_11 = E_k0./(4.*pi.*k0.^4).*(k0.^2 - k1^2 - 2*k1*k30.*xhsi1 + (k1^2 + k2.^2).*xhsi1_q);
end

function phi_22 = phi_22_new(k1,k2,k3,k_vec,beta_vec)
k = sqrt(k1^2 + k2.^2 + k3.^2);

log_beta_vec = interp1(log(k_vec),log(beta_vec),log(k(:)),'linear','extrap');
log_beta = reshape(log_beta_vec,size(k));
beta = exp(log_beta);

k30 = k3 + beta*k1;
k0 = sqrt(k1^2 + k2.^2 + k30.^2);
E_k0 = 1.453*k0.^4./((1 + k0.^2).^(17/6));
C1 = (beta*k1^2).*(k1^2 + k2.^2 - k3.*k30)./(k.^2.*(k1^2 + k2.^2));
C2 = k2.*k0.^2./((k1^2 + k2.^2).^(3/2)).*atan2((beta*k1.*sqrt(k1^2 + k2.^2)),(k0.^2 - k30.*k1.*beta));
xhsi2 = (k2/k1).*C1 + C2;
xhsi2_q = xhsi2.^2;
phi_22 = E_k0./(4.*pi.*k0.^4).*(k0.^2 - k2.^2 - 2.*k2.*k30.*xhsi2 + (k1^2 + k2.^2).*xhsi2_q);
end

function phi_33 = phi_33_new(k1,k2,k3,k_vec,beta_vec)
k = sqrt(k1^2+k2.^2+k3.^2);

log_beta_vec = interp1(log(k_vec),log(beta_vec),log(k(:)),'linear','extrap');
log_beta = reshape(log_beta_vec,size(k));
beta = exp(log_beta);

k30 = k3 + beta*k1;
k0 = sqrt(k1^2+k2.^2+k30.^2);
E_k0 = (1.453*k0.^4./((1+k0.^2).^(17/6)));
phi_33 = (E_k0./(4*pi*(k.^4))).*(k1^2+k2.^2);
end

この図は、元の結果と非常によく一致しているようです。それでも多少の違いはありますが。

補足：シミュレーションではk値100が上限として設定されているため、図のこれより大きい値は正しくありません。それらは、完全な（k2、k3）-「円」のすべての値を使用せずに計算されます。これらの値の偏差も確認できます。

新しい両対数グラフ-F11、F_22、およびF_33のプロット。

于 2013-01-17T09:43:20.070 に答える
1

わかりました、これは私が現時点で得たものです。あなたがそれについてどう思うか聞きたいです-それはまだ完璧ではありません。私は関数にアクセスできないintegralので、（たとえばintegral2私の代わりに）関数を再挿入してコードをテストできれば、現在よりも良い結果が得られる可能性があります。quad2d

私の最初の考えは、-関数のトリプレットごとにforループで計算するbetaことでした。これは非常に遅いことが判明したため、代わりに-valuesのベクトルを使用し、以前と同じように対応するベクトルを計算しました。次に、これら2つのベクトルが渡され、補間関数（）で値が使用されて、特定の値の値が検出されます。[k1,k2,k3]phikbetaphiinterp1betak

function myFunction()

[k1,k2,k3] = deal(0.01:0.1:400);

k = sqrt(k1.^2 + k2.^2 + k3.^2);

c = 1.476;
gamma = 3.9;

Ek = (1.453*k.^4)./((1 + k.^2).^(17/6));
E_int = 1.5-cumtrapz(k,Ek);
beta = c*gamma./(k.*sqrt(E_int));

[F_11,F_22,F_33] = deal(zeros(1,numel(k1)));

k_vec = k;
beta_vec = beta;

for ii = 1:numel(k1)
    phi_11 = @(k2,k3) phi_11_new(k1(ii),k2,k3,k_vec,beta_vec);
    phi_22 = @(k2,k3) phi_22_new(k1(ii),k2,k3,k_vec,beta_vec);
    phi_33 = @(k2,k3) phi_33_new(k1(ii),k2,k3,k_vec,beta_vec);
    F_11(ii) = quad2d(phi_11,-100,100,-100,100);
    F_22(ii) = quad2d(phi_22,-100,100,-100,100);
    F_33(ii) = quad2d(phi_33,-100,100,-100,100);
end

figure
loglog(k1,F_11,'b')
hold on
loglog(k1,F_22,'cyan')
loglog(k1,F_33,'magenta')
hold off
grid on

end

function phi_11 = phi_11_new(k1,k2,k3,k_vec,beta_vec)
k = sqrt(k1^2 + k2.^2 + k3.^2);

beta = reshape(interp1(k_vec,beta_vec,k(:)),size(k));

k30 = k3 + beta*k1;
k0 = sqrt(k1^2 + k2.^2 + k30.^2);
E_k0 = 1.453*k0.^4./((1 + k0.^2).^(17/6));
C1 = (beta*k1^2).*(k1^2 + k2.^2 - k3.*k30)./(k.^2.*(k1^2 + k2.^2));
C2 = k2.*k0.^2./((k1^2 + k2.^2).^(3/2)).*atan2((beta*k1.*sqrt(k1^2 + k2.^2)),(k0.^2 - k30*k1.*beta));
xhsi1 = C1 - (k2/k1).*C2;
xhsi1_q = xhsi1.^2;
phi_11 = E_k0./(4.*pi.*k0.^4).*(k0.^2 - k1^2 - 2*k1*k30.*xhsi1 + (k1^2 + k2.^2).*xhsi1_q);
end

function phi_22 = phi_22_new(k1,k2,k3,k_vec,beta_vec)
k = sqrt(k1^2 + k2.^2 + k3.^2);

beta = reshape(interp1(k_vec,beta_vec,k(:)),size(k));

k30 = k3 + beta*k1;
k0 = sqrt(k1^2 + k2.^2 + k30.^2);
E_k0 = 1.453*k0.^4./((1 + k0.^2).^(17/6));
C1 = (beta*k1^2).*(k1^2 + k2.^2 - k3.*k30)./(k.^2.*(k1^2 + k2.^2));
C2 = k2.*k0.^2./((k1^2 + k2.^2).^(3/2)).*atan2((beta*k1.*sqrt(k1^2 + k2.^2)),(k0.^2 - k30.*k1.*beta));
xhsi2 = (k2/k1).*C1 + C2;
xhsi2_q = xhsi2.^2;
phi_22 = E_k0./(4.*pi.*k0.^4).*(k0.^2 - k2.^2 - 2.*k2.*k30.*xhsi2 + (k1^2 + k2.^2).*xhsi2_q);
end

function phi_33 = phi_33_new(k1,k2,k3,k_vec,beta_vec)
k = sqrt(k1^2+k2.^2+k3.^2);

beta = reshape(interp1(k_vec,beta_vec,k(:)),size(k));

k30 = k3 + beta*k1;
k0 = sqrt(k1^2+k2.^2+k30.^2);
E_k0 = (1.453*k0.^4./((1+k0.^2).^(17/6)));
phi_33 = (E_k0./(4*pi*(k.^4))).*(k1^2+k2.^2);
end

これにより、次の図が生成されます。の最小値では積分は成功しないことに注意してくださいk1。

両対数プロットのF_11、F_22、およびF_33。

編集-ファイ関数内のベータの計算に関するコメント

基本的に私が最初に行ったのと同じことを試したので、-関数beta内でマトリックスを計算する方法の例を追加しましたphi。このコードは非常に遅いため、実際に実行して完了することはありません。

function phi_11 = phi_11_new(k1,k2,k3)
k = sqrt(k1^2 + k2.^2 + k3.^2);

c = 1.476;
gamma = 3.9;
beta = zeros(size(k));
E = @(x) (1.453*x.^4)./((1 + x.^2).^(17/6));
for ii = 1:size(k,1)
    for jj = 1:size(k,2)
        E_int = 1.5-quad(E,0.001,k(ii,jj));
        beta(ii,jj) = c*gamma/(k(ii,jj)*sqrt(E_int));
    end
end


k30 = k3 + beta*k1;
k0 = sqrt(k1^2 + k2.^2 + k30.^2);
E_k0 = 1.453*k0.^4./((1 + k0.^2).^(17/6));
C1 = (beta*k1^2).*(k1^2 + k2.^2 - k3.*k30)./(k.^2.*(k1^2 + k2.^2));
C2 = k2.*k0.^2./((k1^2 + k2.^2).^(3/2)).*atan2((beta*k1.*sqrt(k1^2 + k2.^2)),(k0.^2 - k30*k1.*beta));
xhsi1 = C1 - (k2/k1).*C2;
xhsi1_q = xhsi1.^2;
phi_11 = E_k0./(4.*pi.*k0.^4).*(k0.^2 - k1^2 - 2*k1*k30.*xhsi1 + (k1^2 + k2.^2).*xhsi1_q);
end
于 2013-01-16T14:01:48.873 に答える
0

現在、上記のコードの別のバージョンを実行しています。

それは次のように続きます

clc, clear all, close all
tic



k1 = (0.01:0.1:100);

c = 1.476;
gamma = 3.9;


F_11 = zeros(1,numel(k1));
F_22 = zeros(1,numel(k1));
F_33 = zeros(1,numel(k1));

count = 0;
for i = 1:numel(k1)
    count = count + 1;
    phi_11 = @(k2,k3) phi_11_new(k1,k2,k3,gamma,i);
    phi_22 = @(k2,k3) phi_22_new(k1,k2,k3,gamma,i);
    phi_33 = @(k2,k3) phi_33_new(k1,k2,k3,gamma,i);
    F_11(count) = integral2(phi_11,-100,100,-100,100);
    F_22(count) = integral2(phi_22,-100,100,-100,100);
    F_33(count) = integral2(phi_33,-100,100,-100,100);
end

今回はphi_11、phi_22、phi_33は次のように与えられます

function phi_11 = phi_11_new(k1,k2,k3,gamma,i)
k = sqrt(k1(i).^2 + k2.^2 + k3.^2);
beta = gamma./((k.^(2/3)).*sqrt(hypergeom([1/3,17/6],4/3,-k.^(-2))));  
k30 = k3 + beta.*k1(i);
k0 = sqrt(k1(i).^2 + k2.^2 + k30.^2);
E_k0 = 1.453.*k0.^4./((1 + k0.^2).^(17/6));
C1 = (beta(i).*k1(i).^2).*(k0.^2 - 2.*k30.^2 + beta.*k30.*k1(i))./(k.^2.*(k1(i).^2 + k2.^2));
C2 = k2.*k0.^2./((k1(i).^2 + k2.^2).^(3/2)).*atan2((beta.*k1(i).*sqrt(k1(i).^2 + k2.^2)),(k0.^2 - k30.*k1(i).*beta));
xhsi1 = C1 - k2./k1(i).*C2;
xhsi1_q = xhsi1.^2;
phi_11 = E_k0./(4.*pi.*k0.^4).*(k0.^2 - k1(i).^2 - 2.*k1(i).*k30.*xhsi1 + (k1(i).^2 + k2.^2).*xhsi1_q);
end

function phi_22 = phi_22_new(k1,k2,k3,gamma,i)
k = sqrt(k1(i).^2 + k2.^2 + k3.^2);
beta = gamma./((k.^(2/3)).*sqrt(hypergeom([1/3,17/6],4/3,-k.^(-2)))); 
k30 = k3 + beta.*k1(i);
k0 = sqrt(k1(i).^2 + k2.^2 + k30.^2);
E_k0 = 1.453.*k0.^4./((1 + k0.^2).^(17/6));
C1 = (beta.*k1(i).^2).*(k0.^2 - 2.*k30.^2 + beta.*k1(i).*k30)./(k.^2.*(k1(i).^2 + k2.^2));
C2 = k2.*k0.^2./((k1(i).^2 + k2.^2).^(3/2)).*atan2((beta.*k1(i).*sqrt(k1(i).^2 + k2.^2)),(k0.^2 - k30.*k1(i).*beta));
xhsi2 = (k2./k1(i)).*C1 + C2;
xhsi2_q = xhsi2.^2;
phi_22 = E_k0./(4.*pi.*k0.^4).*(k0.^2 - k2.^2 - 2.*k2.*k30.*xhsi2 + (k1(i).^2 + k2.^2).*xhsi2_q);
end

function phi_33 = phi_33_new(k1,k2,k3,gamma,i)
k = sqrt(k1(i).^2+k2.^2+k3.^2);
beta = gamma./((k.^(2/3)).*sqrt(hypergeom([1/3,17/6],4/3,-k.^(-2)))); 
k30 = k3 + beta.*k1(i);
k0 = sqrt(k1(i).^2+k2.^2+k30.^2);
E_k0 = (1.453.*k0.^4./((1+k0.^2).^(17/6)));
phi_33 = (E_k0./(4*pi.*(k.^4))).*(k1(i).^2+k2.^2);
end

現在、ベータはphi関数内で計算されていることに注意してください。さらに、私はベータの同等の表現を使用しています。詳細については、下の写真をチェックしてください

ここに画像の説明を入力してください

したがって、更新されたモデルでは、パフォーマンスがかなり遅いhypergeomと呼んでいます。これが、最初のコードのように、エネルギースペクトル積分を使用してベータを計算したい主な理由です。

ところで、現時点では、これをうまく実行する方法がわかりません。

よろしく、fpe

于 2013-01-16T10:19:04.317 に答える