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カウンター
カテゴリ
ユーザータグ

LTspiceAkaiKKRmachikaneyamaScilabKKRPSoC強磁性OPアンプPICCPA常微分方程式モンテカルロ解析ecaljodeトランジスタ状態密度インターフェースDOS定電流スイッチング回路PDS5022半導体シェルスクリプトレベルシフト乱数HP6632AR6452AI2C可変抵抗分散関係トランジスタ技術ブレッドボード温度解析反強磁性確率論バンドギャップセミナー数値積分熱設計非線形方程式ソルババンド構造絶縁偏微分方程式ISO-I2CLM358フォトカプラ三端子レギュレータカオスLEDシュミットトリガGW近似A/Dコンバータ発振回路PC817C直流動作点解析USBマフィンティン半径数値微分アナログスイッチTL43174HC4053カレントミラーサーボ量子力学単振り子チョッパアンプ補間2ちゃんねる開発環境bzqltyFFT電子負荷LDAイジング模型BSch基本並進ベクトルブラべ格子パラメトリック解析標準ロジックアセンブラ繰り返し六方最密充填構造SMPコバルトewidthFET仮想結晶近似QSGW不規則合金VCAMaximaGGA熱伝導cygwinスレーターポーリング曲線キュリー温度スイッチト・キャパシタ失敗談ランダムウォークgfortran抵抗相対論位相図スピン軌道相互作用VESTA状態方程式TLP621ラプラス方程式TLP552条件分岐NE555LM555TLP521マントル詰め回路MCUテスタFXA-7020ZR三角波過渡解析ガイガー管自動計測QNAPUPSWriter509ダイヤモンドデータロガー格子比熱熱力学awkブラウン運動起電力スーパーセル差し込みグラフ第一原理計算フェルミ面fsolveCIFxcrysden最大値最小値ubuntu最適化平均場近似OpenMPシュレディンガー方程式固有値問題井戸型ポテンシャル2SC1815TeX結晶磁気異方性OPA2277非線型方程式ソルバフラクタルFSM固定スピンモーメントc/agnuplotPGA全エネルギーfccマンデルブロ集合縮退正規分布キーボード初期値interp1multiplotフィルタ面心立方構造ウィグナーザイツ胞L10構造半金属二相共存ZnOウルツ鉱構造BaOSIC重積分磁気モーメント電荷密度化学反応クーロン散乱岩塩構造CapSenseノコギリ波デバイ模型ハーフメタルRealforceフォノンquantumESPRESSOルチル構造スワップ領域リジッドバンド模型edelt合金等高線凡例軸ラベル線種シンボルトラックボールグラフの分割MAS830LPIC16F785トランス入出力CK1026PC直流解析パラメータ・モデル等価回路モデル不規則局所モーメント関数フィッティング日本語ヒストグラムTS-112ExcelGimp円周率TS-110LMC662片対数グラフ三次元specx.fifortUbuntu文字列疎行列不純物問題ジバニャン方程式ヒストグラム確率論マテリアルデザインP-10境界条件連立一次方程式AACircuit熱拡散方程式HiLAPW両対数グラフ陰解法MBEナイキスト線図負帰還安定性Crank-Nicolson法EAGLE最小二乗法

最新コメント
リンク

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