MATLABで記事の図形を複製する

Question

articleから図を複製したいと思います。

% implementing equation 9 and figure 4 
step = 0.01; t = 1:step:3600; 
d  = 3;  % dimension 
N  = 8000; % number of molecules 
H  = 0.01; % H = [0.01,0.1,1] is in mol/micrometer^3 
H  = H*6.02214078^5; % hence I scaled the Avogadro's number (right or wrong?) 
D  = 10; % diffusion coefficient in micrometer^2/sec 

u(1) = 1./(1.^(d/2)); % inner function in equation 9; first pulse 

for i = 2:numel(t)/1000 
    u(i)  = u(i-1)+(1./(i.^(d/2))); % u-> the pulse number 
    lmda(i) = (1/(4*pi*D))*((N/(H)).*sum(u)).^(2/d); 
end 

figure;plot(lmda) 

しかし、私はそれを複製することはできませんよ。具体的には、私は式9

の表現であるこれまでのところ、私はこのコードが出ていると考えている。図4番を、複製します。パラメータの詳細については

式9

、上記のコードを参照します。著者らは、方程式9における総和はレイマンゼータ系列であることに言及した。それが結果と何か関係があるのだろうか？私は複製しようとしています

図4、：

誰かが親切に私に私が作っています間違いを教えてもらえますか？

P.s：これは宿題ではありません。

Answer 1

問題1：あなたはこの行のアボガドロ数実際に

H = H*6.02214078^5; 

によってスケーリングされ、あなたは約7920=6.022^5によってスケーリングしていると思います。あなたはアボガドロ数でスケールしたいなら、あなたは実行する必要があります。

H = H * 6.02214078e23 % = 6.02214078 * 10^23 : the Avogadro number 

問題2：あなたは本当に（あなたがない限り意味がありません。サンプル数に対してプロットしている、tに対してプロットされていませんがtは整数秒で発生しました）。私たちが何かを見ることができます。この段階では、あなたのループ

for i = 2:numel(t) 
    % ... 
end 
% Then plot 
plot(t, lmda) 

---

から/1000を削除することは本当に間違っています。正しいAvo番号でスケーリングしているので、大規模な注文は途中退出しています。図4のHと式9のHが同じHであることを信頼することをお勧めします。

それに基づいて、間違ったD、N、またはパルス間の時間を使用することをお勧めします。私は以下の私のコードで少し明確にパルスタイミングを設定しました。また、ベクトル化を使用してループをいくらか効率化し、Hスケーリングを削除しました。

dtPulses=100とD=100のように調整すると、プロットはほぼ同じになります。あなたはおそらく、私はパルスが `u`は、一定時間後に生成していないすべての0されているものと

% implementing equation 9 and figure 4 
d = 3;  % dimension 
N = 8000; % number of molecules 
D = 100; % diffusion coefficient in micrometer^2/sec 

dtPulses = 10; % Seconds between pulses 
tPulses = 1:dtPulses:3600; % Time array to plot against 
nt = numel(tPulses); 
i = 1:nt; % pulse numbers 
u = 1 ./ (i.^(d/2)); % inner function in equation 9: individual pulse 
for k = 2:nt % Running sum 
    u(k) = u(k-1)+u(k); 
end 
% Now plot for different H (mol/micrometer^3) 
H = [0.01, 0.1, 1]; 
figure; hold on; linestyles = {':k', '--k', '-k'}; 
for nH = 1:3 
    lmda = ((1/(4*pi*D))*(N/H(nH)).*u).^(2/d);  
    plot(tPulses, lmda, linestyles{nH}, 'linewidth', 2)  
end