差分法による２階常微分方程式の境界値問題の解法：ねむねこ幻想郷：SSブログ

	ブログをはじめるログイン

［境界要素法プログラムを設計する（[計算..｜誤差の蓄積ブログトップ

差分法による２階常微分方程式の境界値問題の解法　[ネコ騙し数学] [編集]

差分法による２階常微分方程式の境界値問題の解法

次の２階常微分方程式

を、差分法を用いて解く解法について考える。

閉区間[a,b]（a≦x≦b）をn等分し

とし、

とおき、さらに分点におけるyの値を

とおく。

（１）式の導関数を

と差分法を用いて近似すると、k=1,2,・・・,n−1において

というn−1本の方程式が得られる。

未知数は、のn−1個なので、連立方程式（４）を解くことにより、これらの未知数を求めることができる。

連立方程式（４）は

とおくと

と書き換えることができる。

そして、（６）式の左辺の係数行列は、３重対角行列なので、TDMA（Tri Diagonal Matrix Algorithm）を用いて解くことができる。

以下に、微分方程式

を解かせたプログラムと、分割数をn=10、n=100とした場合の計算結果を示す。

! 差分法を用いて y''+f(x)y'+g(x)y=h(x) の境界値問題（でぃれくれ型）を解くプログラム
! ただし、固有値問題は解けない!!
parameter (ns=100,ns1=101)
real x(0:ns), y(0:ns)
data x,y/ns1*0.,ns1*0/

n=10

a=0.; b=1. ! 境界のx座標
x(0)=a; x(n)=b
y(0)=2.; y(n)=0. ! 境界条件

call Solver(x,y,n)

write(6,*) ' i      x         y'
do i=0, n
    write(6,100) i, x(i), y(i)
end do

100 format(i3,1x,f9.5,1x,f9.5)
end

function f(x)
f=-14./x
end

function g(x)
g=x**3
end

function h(x)
h=2*x**3
end

! ここより下はいじると危険
! ブラックボックスとして使うべし

subroutine Solver(x,y,n)
real a(n),b(n),c(n),d(n)
real x(0:n), y(0:n)

n1=n-1
dx=(x(n)-x(0))/n

do i=1,n1
    x(i)=x(0)+i*dx
end do

! 差分法によって得られる連立方程式の係数の計算
do i=1,n1
    a(i)=1-dx/2.*f(x(i))
    b(i)=-(2-dx*dx*g(x(i)))
    c(i)=1+dx/2.*f(x(i))
    d(i)=dx*dx*h(x(i))
end do
    d(1)=d(1)-a(1)*y(0) ! 境界条件
    d(n1)=d(n1)-c(n1)*y(n) ! 境界条件

call tdma(a,b,c,d,n1) ! TDMAで連立方程式を解く

do i=1,n1
    y(i)=d(i) ! 計算結果をセット
end do

end

! TDMA
subroutine tdma(a,b,c,d,n)
real a(n), b(n), c(n),d(n)
do i=1,n-1
    ratio=a(i+1)/b(i)
    b(i+1)=b(i+1)-ratio*c(i)
    d(i+1)=d(i+1)-ratio*d(i)
end do
d(n)=d(n)/b(n)
do i=n-1,1,-1
    d(i)=(d(i)-c(i)*d(i+1))/b(i)
end do

end