スポンサーサイト

上記の広告は1ヶ月以上更新のないブログに表示されています。
新しい記事を書く事で広告が消せます。

Scilabでラグランジェ補間 その2

Scilabでラグランジェ補間 その1ではforループを使って多項式補間を行うスクリプトを書きました。しかしながら、Scilabではループ計算は遅いので,できるだけ使わないのがコツなので、今回は多項式生成コマンドであるpolyを利用して高速化を行いました。(LagrangePoly_sce.txt)

clear;

// *** 近似する元の関数 ***
function y = f(x)
// y = 2 .* exp(x - 1) - 1
y = 2 ./ (1 + 9 .* x .^ 2) - 1
endfunction

// *** データ点 ***
N = 10; // データ数
Xn = linspace(-1, 1, N + 1); // xのデータ点
Fn = f(Xn); // yのデータ点
X = linspace(-1, 1, 100); // プロット用のx

// *** 多項式(ラグランジュ)補間 ***
Wn = poly(0,"x");
Pn = 0;
for i = 1:length(Xn) do
Wn = poly(Xn(find(Xn <> Xn(i))),"x");
Pn = Pn + Fn(i) * Wn / horner(Wn,Xn(i));
end

PN = horner(Pn, X);

// *** グラフのプロット ***
subplot(2,1,1);
plot(X, f(X));
plot(Xn, Fn, 'sr');
plot(X, PN, '--k');

// *** 誤差のプロット ***
subplot(2,1,2);
plot(X, PN - f(X));



Scilabでループは遅い


Scilabでラグランジェ補間 その1では任意の関数に対する多項式補間を行う方法を書きました。その際に、多項式を計算するスクリプトとしてScilabで学ぶわかりやすい数値計算法サポートページにて公開されているlagrange.sceを参考にして、以下の様なコードを書きました。

clear;

// *** 近似する元の関数 ***
function y = f(x)
// y = 2 .* exp(x - 1) - 1
y = 2 ./ (1 + 9 .* x .^ 2) - 1
endfunction

// *** データ点 ***
N = 10; // データ数
Xn = linspace(-1, 1, N + 1); // xのデータ点
Fn = f(Xn); // yのデータ点
X = linspace(-1, 1, 100); // プロット用のx

// *** 多項式(ラグランジュ)補間 ***
Pn = zeros(X);
for m = 1:length(X)
for i = 1:(N + 1)
Z(i) = 1;
for j = 1:(N + 1)
if i <> j
Z(i) = Z(i) * (X(m) - Xn(j)) / (Xn(i) - Xn(j));
end
end
Pn(m) = Pn(m) + Fn(i) * Z(i);
end
end

// *** グラフのプロット ***
subplot(2,1,1);
plot(X, f(X));
plot(X, Pn, '--k');
plot(Xn, Fn, 'sr');

// *** 誤差のプロット ***
subplot(2,1,2);
plot(X, Pn - f(X));


しかしながら、このコードはforが3つも入れ子になっています。
通常のプログラミング言語ではあまり問題にならないのかもしれませんが、Scilabに関してはforを多用すると計算速度が遅くなってしまいます。(参考:ループ計算は遅いので,できるだけ使わないScilabでは速いとされる)

そこで本ブログでもScilabで繰り返し計算: ロジスティック写像Scilabの論理演算で条件分岐などループを減らす工夫を行っています。

今回は、多項式の生成コマンドであるpolyを積極的に使って高速化を行いました。

多項式補間(ラグランジェ補間)


N+1 個のデータ点 fn(n=1, 2, …, N+1) を全て通るような多項式は、以下の様に得られます。

数値計算の常識より

\omega (x) = (x - x_1)(x - x_2) \cdots (x - x_{N+1})

\omega_n (x) = \frac{\omega (x)}{x - x_n}

とおくと、求める多項式 PNは以下の様に求められます。

P_N (x) = \sum_{n=1}^{N+1} \frac{\omega_n (x)}{\omega_n (x_n)} f_n

Scilabの数式処理


Scilabは"数値計算"を行うソフトウエアなので、数式を数式の形のままで微分したり、数式同士で割り算をしたり、因数分解したりということは基本的にはやりません。これらの仕事は本来Maximaなどの数式処理ソフトウエアの役目です。

しかしながら、全くそういった仕事ができないかと言うとそういうわけでもなく、多項式に対してこれらの処理を行うことが出来るコマンドを持っています。

今回利用したpolyは、多項式の解から多項式の係数を計算するコマンドです。(参考: scilab 根から多項式を作る!polyscilab 多項式の係数だけを取り出す)

逆に多項式の根を計算するにはrootsを使います。これはScilabで二酸化炭素の状態方程式 その2で利用しました。

多項式を「計算してしまって」数値にするのにはhornerを利用します。

関連エントリ




参考URL




付録


このエントリで使用したScilabのスクリプトを添付します。ファイル名末尾の".txt"を削除して、"_"を"."に変更すれば使えるはずです。(参考:ねがてぃぶろぐの付録)


参考文献/使用機器




フィードバック



にほんブログ村 その他趣味ブログ 電子工作へ

 ↑ 電子工作ブログランキング参加中です。1クリックお願いします。


コメント・トラックバックも歓迎です。 ↓      


 ↓ この記事が面白かった方は「拍手」をお願いします。
スポンサーサイト

tag: Scilab 補間 

Scilabでラグランジェ補間 その1

Scilabでデータの補間ではScilabの標準のコマンドであるinterp1を用いて線形補間やスプライン補間などを行いました。
今回は、補間の考え方の勉強をするために、(実際にはあまり使われない)ラグランジェ補間(多項式補間)をやってみました。

001_20140629231048c59.png

Fig.1: 関数f(x)=2exp(x-1)-1 (青実線)とその3点(赤四角)を通るように計算した多項式補間(黒破線)、及び関数と補間値の誤差(下パネル)。



補間(内挿)とは


数値計算の常識には、内挿に関して以下の様に書かれています。(ただしnの番号付けだけはScilabでスクリプトを書くことを考慮し0からではなく1からに変更しました。これ以降も同様です。)

関数f(x)の値が, 変数xのとびとびの異なる値 x1, x2, …, xN+1に対して与えられているとき, すなわち

fn ≡ f(xn) (n = 1,2, …, N+1)

のみが与えられているとき, xn以外のxに対するf(x)の値を"推測"しようというのが広い意味での"補間"(interpolation; "内挿"ともいう)である.


補間はいろいろな場面で利用されます。例えばScilabでデータの補間AkaiKKRでコバルトの格子定数などです。

今回も例によって数値計算の常識を教科書に補間について勉強します。数値計算の常識の補間に関する大雑把な内容は、以下の通りです。

  1. 最初に考え付くのは多項式(ラグランジェ)補間。しかし、上手く行かない場合も多い。
  2. 次に思いつくのが線形補間。でも折れ曲がりが気になる。
  3. 結局良く使われているのがスプライン補間。


2と3の線形補間とスプライン補間に関してはScilabに標準の命令がありScilabでデータの補間で紹介しています。
そこで今回は(あまり使われない)ラグランジェ補間について書きます。

多項式補間(ラグランジェ補間)


多項式補間の基本的な考え方は N+1 個のデータ点 fn(n=1, 2, …, N+1) を全て通るようなN次多項式は必ず一意に決まるはずなので、それを近似式としましょうと言うものです。

数値計算の常識より

\omega (x) = (x - x_1)(x - x_2) \cdots (x - x_{N+1})

\omega_n (x) = \frac{\omega (x)}{x - x_n}

とおくと、求める多項式 PNは以下の様に求められます。

P_N (x) = \sum_{n=1}^{N+1} \frac{\omega_n (x)}{\omega_n (x_n)} f_n

Scilabスクリプト


Scilabで学ぶわかりやすい数値計算法サポートページにて公開されているlagrange.sceを流用して以下の様なスクリプトを作成しました。(Lagrange_sce.txt)
clear;

// *** 近似する元の関数 ***
function y = f(x)
// y = 2 .* exp(x - 1) - 1
y = 2 ./ (1 + 9 .* x .^ 2) - 1
endfunction

// *** データ点 ***
N = 10; // データ数
Xn = linspace(-1, 1, N + 1); // xのデータ点
Fn = f(Xn); // yのデータ点
X = linspace(-1, 1, 100); // プロット用のx

// *** 多項式(ラグランジュ)補間 ***
Pn = zeros(X);
for m = 1:length(X)
for i = 1:(N + 1)
Z(i) = 1;
for j = 1:(N + 1)
if i <> j
Z(i) = Z(i) * (X(m) - Xn(j)) / (Xn(i) - Xn(j));
end
end
Pn(m) = Pn(m) + Fn(i) * Z(i);
end
end

// *** グラフのプロット ***
subplot(2,1,1);
plot(X, f(X));
plot(X, Pn, '--k');
plot(Xn, Fn, 'sr');

// *** 誤差のプロット ***
subplot(2,1,2);
plot(X, Pn - f(X));


結果


近似する関数は数値計算の常識に倣って以下の2種類を行いました。

f(x)=2 \exp(x-1) - 1

f(x)=\frac{2}{1+9x^2}-1

定義域はともに -1 ≦ x ≦ 1 です。
以下に示すように、前者は上手く補間できるのですが、後者はイマイチです。

002_20140629231048cb2.png

Fig.2: 関数f(x)=2exp(x-1)-1 (青実線)とその7点(赤四角)を通るように計算した多項式補間(黒破線)、及び関数と補間値の誤差(下パネル)。Fig.1と比較して、通るべき点の数を増やすことによって近似の誤差が少なくなっていることが読み取れる。


003_20140629231048721.png

Fig.3: 関数 2/(1+9x2)-1 (青実線)とその11点(赤四角)を通るように計算した多項式補間(黒破線)、及び関数と補間値の誤差(下パネル)。通るべき点の数を増やしても近似の誤差が小さくならない問題点がある。


Fig.1-2は補間が上手く行っている例です。しかしながら、Fig.3の様に補間するもとの関数によっては誤差が大きくなってしまうことがあるのがラグランジェ補間の問題点です。

スプライン補間


もう少し優等生なスプライン補間を用いた結果も比較のため載せておきます。(spline_sce.txt)

004_20140629233541671.png

Fig.4: スプライン補間。図の読み方はFig.3と同じ。Fig.3と比較して誤差が少なくなっていることが分かる。


Fig.3とFig.4を比較するとスプライン補間のほうが良い結果を示していることがわかります。
ただ、「スプライン補間さえ使ってればいつでもオッケー」と言うわけではないと言うことも覚えておかなければいけません。

数値計算の常識曰く

どだい, f(xn) (n = 1, 2, …, N+1) から f(x) (x≠xn)を知ろうというのは, 関数fがよっぽど良い性質を持っていなければ"原理的に"無理であるのは明らかであろう. (x = x1, …, xN+1 で一致するけれど他のxでは値が異なる関数はいくらでもありうる.)


関連エントリ




参考URL




付録


このエントリで使用したScilabのシミュレーション用ファイルを添付します。ファイル名末尾の".txt"を削除して、"_"を"."に変更すれば使えるはずです。(参考:ねがてぃぶろぐの付録)


参考文献/使用機器




フィードバック



にほんブログ村 その他趣味ブログ 電子工作へ

 ↑ 電子工作ブログランキング参加中です。1クリックお願いします。


コメント・トラックバックも歓迎です。 ↓      


 ↓ この記事が面白かった方は「拍手」をお願いします。

tag: Scilab 補間 

AkaiKKRでNiCuの磁気モーメント

AkaiKKRでFeCoの磁気モーメントと格子定数AkaiKKRでFeCrの磁気モーメントに続いて不規則ニッケル銅合金の第一原理計算をAkaiKKR(Machikaneyama)を利用して計算しました。

NiCu.gif

Fig.1: fcc-Ni100-xCux合金の状態密度の銅濃度依存性。アニメーションは純ニッケルから順に10%刻みで純銅まで。



スレーターポーリング曲線


AkaiKKRでFeCoの磁気モーメントと格子定数AkaiKKRでFeCrの磁気モーメントとスレーターポーリング曲線の中の頂点と左端の計算をAkaiKKR(Machikaneyama)を用いて行ってきました。今回は、このスレーターポーリング曲線の右端に位置するNiCu合金の第一原理計算を行います。

以下に示すのは、磁石・磁気の用語辞典にて紹介されているスレーターポーリング曲線です。

001_20140604225835b78.jpg

Fig.2: スレーターポーリング(Slater-Pauling)曲線


計算手法


AkaiKKRでFeCoの磁気モーメントと格子定数AkaiKKRでFeCrの磁気モーメントと同様に格子定数と化学組成を変化させ、それぞれの組成に対してMurnaghanの状態方程式へ全エネルギーをフィッティングします。フィッティングから得られた平衡格子定数での磁気モーメントは、スプライン関数を用いて内挿しました。

結晶構造は全て面心立方構造(fcc)とし、格子定数の計算範囲は 6.50 ≦ a ≦ 6.80 で 0.02 bohr 刻みで変化させています。入力ファイルの例は以下の通りです。

 go    data/nicu_6.64
fcc 6.64 , , , , , ,
0.001 1.2 sra mjw mag 2nd
update 16 200 0.02
1
NiCu 2 0 0 2 28 40
29 60
1
0.00000 0.00000 0.00000 NiCu


計算結果と議論


得られた磁気モーメントは、スレーターポポーリング曲線を再現しNi40Cu60組成付近で、磁気モーメントを失いました。

001_201406281825354ad.png

Fig.3: 格子定数と磁気モーメントの銅濃度依存性


これに対して状態密度は興味深い結果が得られました。
ニッケルと銅は共にfcc構造であり、周期表でも隣り合わせの関係にあるので、リジッドバンドモデル的な状態密度の変化をするのかと思いきや、銅の濃い組成では状態密度の形が鋭さをずいぶんと失ってしまっています。

Ni100Cu0.png

Fig.4: Ni100Cu0


Ni90Cu10.png

Fig.5: Ni90Cu10


Ni80Cu20.png

Fig.6: Ni80Cu20


Ni70Cu30.png

Fig.7: Ni70Cu30


Ni60Cu40.png

Fig.8: Ni60Cu40


Ni50Cu50.png

Fig.9: Ni50Cu50


Ni40Cu60.png

Fig.10: Ni40Cu60


Ni30Cu70.png

Fig.11: Ni30Cu70


Ni20Cu80.png

Fig.12: Ni20Cu80


Ni10Cu90.png

Fig.13: Ni10Cu90


Ni0Cu100.png

Fig.14: Ni0Cu100


この計算の結果がどの程度現実の銅ニッケル合金を再現しているのかは私には判断できませんが、不規則性の効果が最も強く現れるのは比較的深いエネルギー準位の部分で、現実の物性の大部分を占めると思われるフェルミ準位近傍では支配的ではないかもしれません。

関連エントリ




参考URL




参考文献/使用機器




フィードバック



にほんブログ村 その他趣味ブログ 電子工作へ

 ↑ 電子工作ブログランキング参加中です。1クリックお願いします。


コメント・トラックバックも歓迎です。 ↓      


 ↓ この記事が面白かった方は「拍手」をお願いします。

tag: AkaiKKR machikaneyama KKR CPA 強磁性 スレーターポーリング曲線 状態密度 DOS Scilab 

AkaiKKRでFeCrの磁気モーメント

AkaiKKRでFeCoの磁気モーメントと格子定数に引き続きAkaiKKR(Machikaneyama)を用いてスレーターポーリング曲線に乗る強磁性遷移金属合金のひとつである不規則FeCr系合金の磁気モーメントを計算しました。

FeCr.gif

Fig.1: bcc-FexCr100-x合金の状態密度のクロム濃度依存性。アニメーションは純鉄から順に10%刻みで純クロムまで。


得られた磁気モーメントは、スレーターポーリング曲線が示すとおり、平均価電子数に対して直線的に変化しました。この結果は、リジッドバンドモデルを用いて理解することが出来ます。ただし、コヒーレントポテンシャル近似を用いて得られた状態密度曲線は、不規則性の効果により、かなり鋭さを減じています。


スレーターポーリング曲線


AkaiKKRでFeCoの磁気モーメントと格子定数ではAkaiKKR(Machikaneyama)を用いてスレーターポーリング(Slater-Pauling)曲線の頂点に位置するFeCo系合金の計算を行いました。
今回は、スレーターポーリング曲線の左端で強磁性が消失する部分に相当するFeCr合金の磁気モーメントを計算します。

以下に示すのは、磁石・磁気の用語辞典にて紹介されているスレーターポーリング曲線です。

001_20140604225835b78.jpg

Fig.2: スレーターポーリング(Slater-Pauling)曲線


計算手法


AkaiKKRでFeCoの磁気モーメントと格子定数と同様に格子定数と化学組成を変化させ、それぞれの化学組成に対してMurnaghanの状態方程式をフィッティングし、平衡格子定数を求めた後、スプライン補間から平衡格子定数のときの磁気モーメントを求めました。

結晶構造は全て体心立方構造(bcc)とし、格子定数の範囲は a = 5.23 bohr から a = 5.35 bohr まで 0.01 bohr 刻みとしました。以下は、FeCr系入力ファイルの一例です。

 go    data/fecr_5.29
bcc 5.29 , , , , , ,
0.001 1.2 sra mjw mag 2nd
update 16 200 0.02
1
FeCr 2 0 0 2 26 30
24 70
1
0.00000 0.00000 0.00000 FeCr


結果と議論


得られた磁気モーメントは、スレーターポーリング曲線を再現し、純Crで強磁性が消失しました。

001_20140605072505cfc.png

Fig.3: 磁気モーメントの鉄濃度依存性


次に状態密度を見てみます。
まずはコヒーレントポテンシャル近似(CPA)が必要ない端成分である純鉄と純クロムの状態密度を以下に示します。

Fe100Cr0.png
Fig.4: 純鉄の状態密度

Fe0Cr100.png
Fig.5: 純クロムの状態密度


クロムは非磁性なので、多数スピンバンドと少数スピンバンドの状態密度が同じ形状をしています。
しかしながら、各スピンバンドごとの状態密度曲線の形状は、鉄とクロムで互いによく似ています。このことからリジッドバンドモデルが成り立つのではないかと言う予測が立ちます。

リジッドバンド的な考え方では、FeCr系の磁気モーメントの変化は、鉄から価電子数を減らしていったとき、多数スピンバンドの状態密度がフェルミ準位に対して、相対的に高エネルギー側へシフトし、少数スピンバンドではほとんど変わらないと解釈されます。

このことに関して、コヒーレントポテンシャル近似による合金化の効果も含めて状態密度を順に示します。

Fe100Cr0.png
Fig.6: Fe100Cr0の状態密度

Fe90Cr10.png
Fig.7: Fe90Cr10の状態密度

Fe80Cr20.png
Fig.8: Fe80Cr20の状態密度

Fe70Cr30.png
Fig.9: Fe70Cr30の状態密度

Fe60Cr40.png
Fig.10: Fe60Cr40の状態密度

Fe50Cr50.png
Fig.11: Fe50Cr50の状態密度

Fe40Cr60.png
Fig.12: Fe40Cr60の状態密度

Fe30Cr70.png
Fig.13: Fe30Cr70の状態密度

Fe20Cr80.png
Fig.14: Fe20Cr80の状態密度

Fe10Cr90.png
Fig.15: Fe10Cr90の状態密度

Fe0Cr100.png
Fig.16: Fe0Cr100の状態密度


CPA計算の結果から、多数スピンバンドの状態密度は、不規則性の効果によって形状の鋭さがかなり失われています。これに対して、少数スピンバンドはあまり変化しているように見えません。
FeCr系の状態密度は常に多数スピンバンドのフェルミ準位がdバンドの中にあり、強磁性となるときは弱い強磁性です。

関連エントリ




参考URL




参考文献/使用機器




フィードバック



にほんブログ村 その他趣味ブログ 電子工作へ

 ↑ 電子工作ブログランキング参加中です。1クリックお願いします。


コメント・トラックバックも歓迎です。 ↓      


 ↓ この記事が面白かった方は「拍手」をお願いします。

tag: AkaiKKR machikaneyama KKR CPA 強磁性 Scilab スレーターポーリング曲線 DOS 状態密度 

FC2カウンター
カテゴリ
ユーザータグ

LTspiceAkaiKKRScilabmachikaneyamaKKRPSoCCPAOPアンプPIC強磁性モンテカルロ解析常微分方程式トランジスタodeインターフェース状態密度DOSecalj定電流PDS5022スイッチング回路半導体シェルスクリプト乱数レベルシフトHP6632A温度解析ブレッドボードI2CR6452A分散関係トランジスタ技術可変抵抗確率論数値積分反強磁性セミナー非線形方程式ソルバ絶縁バンドギャップ熱設計偏微分方程式バンド構造GW近似カオス三端子レギュレータLEDフォトカプラシュミットトリガISO-I2CA/DコンバータLM358USBカレントミラーTL431マフィンティン半径PC817C数値微分アナログスイッチ発振回路サーボ直流動作点解析74HC40532ちゃんねる標準ロジックチョッパアンプLDAアセンブラFFTbzqltyイジング模型ブラべ格子開発環境補間量子力学電子負荷BSchパラメトリック解析単振り子基本並進ベクトル熱伝導繰り返しGGAMaximaTLP621ewidthSMP相対論抵抗位相図ランダムウォークスピン軌道相互作用六方最密充填構造不規則合金FETコバルト失敗談QSGWcygwinスレーターポーリング曲線スイッチト・キャパシタラプラス方程式gfortranキュリー温度状態方程式条件分岐格子比熱TLP552LM555TLP521三角波NE555過渡解析FXA-7020ZRWriter509テスタ詰め回路MCUマントルダイヤモンドQNAPデータロガーガイガー管自動計測UPS井戸型ポテンシャルawk第一原理計算仮想結晶近似ブラウン運動差し込みグラフ平均場近似fsolve起電力熱力学OpenMPスーパーセル固有値問題最適化最小値VCAシュレディンガー方程式VESTAubuntu最大値面心立方構造PGAOPA2277L10構造非線型方程式ソルバ2SC1815fccフェルミ面等高線ジバニャン方程式ヒストグラム確率論マテリアルデザイン正規分布結晶磁気異方性interp1フィルタ初期値ウィグナーザイツ胞c/aルチル構造岩塩構造スワップ領域リジッドバンド模型edeltBaOウルツ鉱構造重積分SIC二相共存ZnOquantumESPRESSOCapSensegnuplotmultiplot全エネルギー固定スピンモーメントFSM合金ノコギリ波フォノンデバイ模型ハーフメタル半金属TeXifortTS-110不規則局所モーメントTS-112等価回路モデルパラメータ・モデルヒストグラムExcel円周率GimpトラックボールPC直流解析入出力文字列マンデルブロ集合キーボードフラクタル化学反応三次元Realforce縮退日本語最小二乗法関数フィッティング疎行列シンボル線種ナイキスト線図陰解法負帰還安定性熱拡散方程式EAGLECrank-Nicolson法連立一次方程式P-10クーロン散乱Ubuntu境界条件MBEHiLAPW軸ラベルトランスCK1026MAS830L凡例PIC16F785LMC662AACircuit両対数グラフ片対数グラフグラフの分割specx.f

最新コメント
リンク

にほんブログ村 その他趣味ブログ 電子工作へ
上記広告は1ヶ月以上更新のないブログに表示されています。新しい記事を書くことで広告を消せます。