Scilabで熱拡散方程式 その3 (陰解法)

Scilabで熱拡散方程式 その1 (陽解法)Scilabで熱拡散方程式 その2 (無次元化)では陽解法で偏微分方程式を解く方法を書くと同時に、時間と空間の刻み幅を自由に選べないという問題点を指摘しました。今回はその問題点を解決した陰解法(Crank-Nicolson法)を紹介します。


001_20140215031531793.png
Fig.1: 一次元熱伝導問題の解



偏微分方程式の陰解法


Scilabで熱拡散方程式 その1 (陽解法)では、Octaveの精義を参考にして、一次元の熱伝導を題材に偏微分方程式を陽解法で解くということを行いました。しかしながらScilabで熱拡散方程式 その2 (無次元化)で紹介したとおり、陽解法では λ ≦ 1/2 という安定性の条件から空間の分解能と時間の分解能を独立に決めることが出来ないという欠点があることが確認されました。

Octaveの精義では、この問題を陰解法のひとつであるCrank-Nicolson法を用いて解決する方法を紹介しています。そこで今回は、その1の計算をCrank-Nicolson法で計算します。

計算


計算結果が冒頭に示したFig.1です。
そのスクリプトはCN_sce.txtです。

clear;

// *** 時間と空間の分割 ***
n = 129; // 空間の分割数
dx = 1 / (n - 1); // 空間の刻み幅

r = 4

dt = r * dx ^ 2; // 時間の刻み幅
m = 1001; // 時間の分割数
dm = 100; // プロットする m の間隔

// *** 初期条件と境界条件 ***
// 初期条件
u = zeros(n,m);
// 境界条件
u(1,1) = 0;
u(n,1) = 1;

// *** 行列C ***
C = sparse(toeplitz([-2, 1, zeros(1,n - 2)]));
A = sparse(2 * eye(n,n) - r * C);
B = sparse(2 * eye(n,n) + r * C);
b = zeros(n,1);
b($) = 2 * r;

// *** 偏微分方程式の計算 ***
for j = 1:(m - 1)
// 境界以外の計算
u(:, j + 1) = A \ (B * u(:,j) + b);
// 境界条件
u(1, j + 1) = 0;
u(n, j + 1) = 1;
end

// *** グラフのプロット ***
plot([0:dx:1],u(:,1:dm:m));


陰解法の利点


Scilabで熱拡散方程式 その1 (陽解法)Scilabで熱拡散方程式 その2 (無次元化)では、陽解法の問題点として時間と空間の刻み幅を独立に決定することが出来ないという点を挙げました。

Crank-Nicolson法はこの問題点を改善したものです。


002_201402150315301ea.png
Fig.2: 1次元熱伝導方程式の陽解法(FTCS法)による数値解で、Octaveの精義のコードをScilabへ移植したもの。全てh2=1/322に固定しr = 1/6 (左上), 1/2 (右上), 65/128 (左下)と変化させた結果。r = 64/128 > 1/2 の場合には解が不安定になっていることが読み取れる。

003_20140215031530c4c.png
Fig.3: 1次元熱伝導方程式の陰解法(Crank-Nicolson法)による数値解で、Octaveの精義のコードをScilabへ移植したもの。全てh2=1/1282に固定しr=1 (左上),2 (右上), 4(左下)と変化させた結果。全てのrに対して安定な解が得られていることが分かる。



関連エントリ




参考URL




付録


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


参考文献/使用機器






フィードバック



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

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


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


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

tag: Scilab 偏微分方程式 熱伝導 陰解法 Crank-Nicolson法 

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

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

最新コメント
リンク

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