So-net無料ブログ作成

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

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


 


次の2階常微分方程式



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


 


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



とし、



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



とおく。


 


(1)式の導関数を



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



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


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


 


連立方程式(4)は



とおくと



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


そして、(6)式の左辺の係数行列は、3重対角行列なので、TDMATri Diagonal Matrix Algorithm)を用いて解くことができる。


 


以下に、微分方程式



を解かせたプログラムと、分割数をn=10n=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



 


計算結果はおかしく見えるかもしれないけれど、この計算結果は正しいケロよ。


 


上の微分方程式の解は、ベッセル関数という特殊関数を用いないと表わせないうえに、ベッセル関数は組み込み関数にないので、n=100のときの数値解を厳密解だと考えて欲しいにゃ。


 


少し違うけれど、十進BASIC用のプログラムは 
http://nekodamashi-math.blog.so-net.ne.jp/2016-01-24


に書いてあるので、Fortranを使えないヒト、あるいは、Fortranから他のコンピュータ言語にできないヒトは、コチラのプログラムを参考にしてほしいにゃ。


 


タグ:数値解析
nice!(0)  コメント(0) 

nice! 0

コメント 0

コメントを書く

お名前:
URL:
コメント:
画像認証:
下の画像に表示されている文字を入力してください。

この広告は前回の更新から一定期間経過したブログに表示されています。更新すると自動で解除されます。