差分法で熱伝導方程式を解く

差分法では、偏導関数を次のように近似する。

　（前進差分）もしくは、
　（後退差分）
2階の偏微分は1階の偏導関数の差分であるから、前進差分と後退差分を用いれば、次式が得られる。

　（中心差分）

---

熱伝導方程式

差分法を用いて熱伝導方程式の数値解を求めてみる。一次元の熱伝導方程式は次式で表される。


Tは温度 [K]、αは温度伝導率（温度拡散率）であり、次式より求まる。

ここで、λは熱伝導率 [W/(m・K)]、ρは密度 [kg/m³]、c_pは定圧比熱 [J/(kg・K)]である。

---

陽解法

前進差分を用いて数値解を求める方法は陽解法と呼ばれる。T_x,tとして、陽解法による差分近似式は次式となる。


T_i,j+1について解けば、


が得られる。ここで、γ=Δt/Δx²である。この式より、T_i-1,j, T_i,j, T_i+1,jの3点から、次の時刻のT_i,j+1を計算することができる。

Pythonでプログラムを書くとこうなる。

#coding: utf-8

L = 0.1  #長さ[m]
M = 10   #分割数
dx = L/M #空間刻み[m]

dt = 1.0 #時間刻み[s]
N = 100  #時間ステップ数

#温度伝導率(鉄)
alpha = 80.2/(7874.0*440.0)

gamma = dt/dx**2
a = alpha*gamma

#初期条件(温度273.15K)
T = [273.15 for i in range(M+1)]

#境界条件
T[0] = 300.0  #温度固定
T[M] = T[M-1] #断熱

f = open('output','w')

for j in range(1,N):
    for i in range(1,M):
        preT = T

        #陽解法による差分式
        T[i] = a*preT[i+1]+(1-2*a)*preT[i]+a*preT[i-1]
        T[0] = 300.0
        T[M] = T[M-1]

    #計算結果をファイルへ出力
    f.write('# t=%ss\n'%(j*dt))
    for i in range(M+1):
        f.write('%s, %s\n'%(i*dx,T[i]))
    f.write('\n\n')

f.close()

グラフにプロットするとこんな感じ。時間とともに、左側から温度が上昇していく様子がわかる。

---

陰解法

後退差分を用いて数値解を求める方法は陰解法と呼ばれる。差分近似式は次式となる。


整理すると、


陽解法とは異なり、T_i,jから、次の時刻でのT_i-1,j+1, T_i,j+1, T_i+1,j+1の3点を求めなければならない。前式は行列で書くと次のようになる。


ただし、行列のT_1,j+1とT_M,j+1には境界条件を設定する。従って、T_{(1,2,...,M),j+1}は次式より求まる。


Python + numarray でプログラムを書くとこうなる。

#coding: utf-8

import numarray
import numarray.linear_algebra as la

L = 0.1  #長さ[m]
M = 10   #分割数
dx = L/M #空間刻み[m]

dt = 1.0 #時間刻み[s]
N = 100  #時間ステップ数

#温度伝導率(鉄)
alpha = 80.2/(7874.0*440.0)

gamma = dt/dx**2
a = alpha*gamma

#初期条件(温度273.15K)
T = numarray.array([[273.15,] for i in range(M+1)])
#境界条件
T[0][0] = 300.0  #温度固定
T[M][0] = T[M-1][0] #断熱

#行列の作成
A = numarray.array([[0. for i in range(M+1)] for j in range(M+1)])
for i in range(1,M):
    A[i][i-1] = -a
    A[i][i]   = 1+2*a
    A[i][i+1] = -a
#境界条件
A[0][0] = 1.   #一定温度
A[M][M]   = 1. #断熱

#Aの逆行列
A_inv = la.inverse(A)

f = open('output','w')

for j in range(1,N):
    preT = T

    T = numarray.dot(A_inv, preT)
    T[0][0] = 300.0
    T[M][0] = T[M-1][0]

    #計算結果をファイルへ出力
    f.write('# t=%ss\n'%(j*dt))
    for i in range(M+1):
        f.write('%s, %s\n'%(i*dx,T[i][0]))
    f.write('\n\n')

f.close()

グラフにプロットするとこんな感じ。

---

クランク・ニコルソン法

ルンゲ・クッタ法

ルンゲ・クッタ法を用いて、微分方程式の初期値問題を解く。

微分方程式と初期条件が次式で与えられたとする。




このとき、4次のルンゲ・クッタ法により次式が与えられる。


ここで h は刻み幅、










である。初期条件x₀，y₀からx₁，y₁が求まり、x₁，y₁からx₂，y₂が求まる・・・というふうに計算を繰り返して数値解を得る。

例えば、微分方程式 f(x,y) = y-3x と、初期条件 y(0)=1 が与えられたとすると、

#coding: utf-8

def f(x, y):
    return (y-3*x)

#初期条件
x=0.
y=1.

#増加幅
h=0.1

while x <= 1:
    print x,y

    k1 = f(x, y)
    k2 = f(x+h/2, y+h/2*k1)
    k3 = f(x+h/2, y+h/2*k2)
    k4 = f(x+h, y+h*k3)
    k  = (k1+2*k2+2*k3+k4)/6
    y = y+k*h
    x+=h

計算すると、

x   y
0.0 1.0
0.1 1.08965833333
0.2 1.1571948583
0.3 1.20028300587
0.4 1.21635151984
0.5 1.20255872281
0.6 1.15576407582
0.7 1.07249674681
0.8 0.948920873415
0.9 0.78079717244
1.0 0.56344051173

と求まる。

ちなみに解析解は y=3x-2e^x+3 である。グラフにプロットするとこんな感じ。良い精度で近似できていることがわかる。

Python - numarray で最小二乗法

Python - numarray　で最小二乗法をやってみる。

当てはめる関数を
とする。

係数a₁,a₂, ... , a_kを求めるには、正規方程式
を解けばよい。

ここで、


　，　である。


例えば、

(x,y) = (1.0, 8.1), (2.0, 2.5), (3.0, 2.1), (4.0, 3.4), (5.0, 9.8)
というデータがあったとして、これに二次関数


を当てはめることを考える。

, ,
である。

Pythonでnumarrayを使って行列演算をする。


# coding: utf-8
import numarray
import numarray.linear_algebra as la

F = numarray.array([[1,1,1],[4,2,1],[9,3,1],[16,4,1],[25,5,1]])
y = numarray.array([[8.1],[2.5],[2.1],[3.4],[9.8]])

tF = numarray.transpose(F)
tFF = numarray.dot(tF, F)
tFF_inv = la.inverse(tFF) #tFFの逆行列
tFy = numarray.dot(tF, y)
a = numarray.dot(tFF_inv, tFy)

print a


計算してみると、(a₁, a₂, a₃) = (1.8357, -10.5843, 16.74)　と求まる。実際にプロットするとこんな感じ。

SQLite - Python

ちょっとしたデータベースを作りたくて、pythonからSQLiteを使ってみる。ライブラリリファレンスが翻訳されてないなぁ。とりあえずサンプル。

#coding: utf-8

import sqlite3

#データベースへ接続（存在していなければ作成）
con = sqlite3.connect("test.db")

#カーソルって何？
cur = con.cursor()

#本のデータを収める'book'テーブルを作成
cur.execute("""
create table if not exists booklist(
id intager not null primary key,
title text,
author text);
""")

#テーブルにデータを挿入
cur.execute("""
insert into booklist(id, title, author)
values(1, 'Python Tutorial', 'Guide van Rossum');
""")

cur.execute("""
insert into booklist(id, title, author)
values(2, 'Microeconomics', 'Paul Krugman, Robin Wells');
""")

#データを取り出す
cur.execute("select * from booklist;")
print cur.fetchall()

cur.execute("select title from booklist;")
print cur.fetchall()

cur.execute("select * from booklist where id=1;")
print cur.fetchall()

con.close()

ニュートン法

学生時代にポケットコンピュータのBASICで組んだやつを見つけて、Pythonで書いてみた。なつかし。

#coding: utf8

print "「ax^2+bx+c=0」の解をニュートン法により求める"
a = float(input('a= '))
b = float(input('b= '))
c = float(input('c= '))

x = 1
n = 0

while 1:
    f = a*x*x+b*x+c
    df= a*x+b
    xn = x-f/df
    if(abs((x-xn)/x)<1e-10):
        print "解はx=",xn
        break
    else: x=xn

    n+=1
    if n>1e5:
        print "解なし"
        break