从这段代码开始:
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
这个过程导致我得到的结果与其他研究得出的结果不匹配。我应该匹配的结果发布如下:
而我的看起来像这样
很容易看出只有比较结果才符合理论结果;因此,我认为该缺陷可能存在于函数 phi_11_new(和 phi_22_new)之外的 beta 定义中。
你们中的任何人都可以建议如何在 phi_11_new(和 phi_22_new)中计算 beta,而不是像我现在那样?
提前感谢大家的支持。
最诚挚的问候, FPE
最佳答案
我改进了插值,使其不再因小值而崩溃。它还返回更正确的值,因为它现在插入值的对数。
这是现在的代码。
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)-“圆”中的所有值。我们还可以看到这些值的偏差。
关于function - 构造一个函数,我们在Stack Overflow上找到一个类似的问题: https://stackoverflow.com/questions/14337208/