偏微分方程式解法

差分法

差分近似式の導出を次にしめします。

2次元非定常解析だけ示します。

はじめに、差分近似式を示します。

そして、問題文を書きます。

その後、Excel macroの計算結果とExcel macroを示します。

Excel　macroのコイピー＆ペイストの方法はDFTと同じです。

Sheet1とSheet2を用意（作る）してください。

実行すれば、計算結果が表示（数値のみ）されます。

後は、ご自身でグラフを作成してください。

Excelは、ここまでです。

以下、FORTRANの説明です。

ここに出ているFORTRANで関数のプログラムをコピー＆ペイストすると、1行おきにブランク行が出来るが、構わず、******.forに保存して、コンパイル、実行出来ることを確認しました。

謝辞

偏微分方程式の導出は、工業数学＜上＞を参考にさせて頂きました。わかりやすい説明で理解の

助けになりました。

C.R.ワイリー著　工業数学　富久泰明　訳

例題　一次元非定常熱伝導問題

一次元の棒を考えてください。

初期温度100℃に全体がなっていました。

時刻ｔ＝0で両端が0℃になりました。

フーリエ数：Fo＝ｃρ/λL＾2

ここで、ｃ：比熱

　　　　ρ：密度

　　　　λ：熱伝導率

　　　　L：棒の長さ

フーリエ数：Fo＝0.00125、0.0125、0.0625、0.125、0.25

の時の温度分布を計算します。

サンプルプログラムをつぎにしめします。

0．Excel macro　FDM

      Sub FDM()

'

'     FDM.xlsm

'

      Dim temp(101, 3) As Single

      Dim tempnew(101, 3) As Single

      Dim resultemp(101, 11) As Single

      Dim Fo(5) As Single

      Dim IPTIME(5) As Integer

      Dim temp1(5) As Single

      Dim dt As Single

      Dim time As Single

      Dim itime As Integer

      dx = 0.001

      dy = 0.001

      dt = 0.01

      SL = dx * 100#

      ramda = 46#

      c = 460.5

      roh = 7870#

      Fo(1) = 0.00125

      Fo(2) = 0.0125

      Fo(3) = 0.0625

      Fo(4) = 0.125

      Fo(5) = 0.25

      AKAPPA = ramda / c / roh

      For i = 1 To 5

      IPTIME(i) = Fo(i) / AKAPPA * SL * SL / dt

      Next i

      A1 = dt * ramda / c / roh / dx / dx

      A2 = A1

      A3 = dt * ramda / c / roh / dy / dy

      A4 = A3

      AA = 1# - (A1 + A2 + A3 + A4)

      Sheets("Sheet1").Select

      ActiveSheet.[A1].Select

      ActiveCell.Offset(0, 0) = "time"

      ActiveCell.Offset(0, 1) = "FDM:x=0.01m"

      ActiveCell.Offset(0, 2) = "FDM:x=0.02m"

      ActiveCell.Offset(0, 3) = "FDM:x=0.03m"

      ActiveCell.Offset(0, 4) = "FDM:x=0.04m"

      ActiveCell.Offset(0, 5) = "FDM:x=0.05m"

      ActiveCell.Offset(0, 6) = "theory:x=0.01m"

      ActiveCell.Offset(0, 7) = "theory:x=0.02m"

      ActiveCell.Offset(0, 8) = "theory:x=0.03m"

      ActiveCell.Offset(0, 9) = "theory:x=0.04m"

      ActiveCell.Offset(0, 10) = "theory:x=0.05m"

      For i = 1 To 101

      For j = 1 To 3

      temp(i, j) = 100#

      Next j

      Next i

      temp(1, 1) = 0#

      temp(101, 1) = 0#

      temp(1, 2) = 0#

      temp(101, 2) = 0#

      temp(1, 3) = 0#

      temp(101, 3) = 0#

      For itime = 1 To IPTIME(5) + 3

      time = itime * dt

      j = 1

      For i = 2 To 100

      tempnew(i, j) = A1 * temp(i + 1, j) + A2 * temp(i - 1, j) + A3 * temp(i, j + 1) + A4 * temp(i, j + 1) + AA * temp(i, j)

      Next i

      j = 2

      For i = 2 To 100

      tempnew(i, j) = A1 * temp(i + 1, j) + A2 * temp(i - 1, j) + A3 * temp(i, j - 1) + A4 * temp(i, j + 1) + AA * temp(i, j)

      Next i

      j = 3

      For i = 2 To 100

      tempnew(i, j) = A1 * temp(i + 1, j) + A2 * temp(i - 1, j) + A3 * temp(i, j - 1) + A4 * temp(i, j - 1) + AA * temp(i, j)

      Next i

      For i = 1 To 5

      If (IPTIME(i) = itime) Then

      For k = 1 To 101

      resultemp(k, i + 1) = tempnew(k, 2)

      Next k

      Else

      End If

      Next i

      For i = 1 To 101

      For j = 1 To 3

      temp(i, j) = tempnew(i, j)

      Next j

      Next i

      Call theory_solution_time(time, temp1, ramda, c, roh, SL)

      ActiveCell.Offset(itime, 0) = time

      ActiveCell.Offset(itime, 1) = temp(11, 2)

      ActiveCell.Offset(itime, 2) = temp(21, 2)

      ActiveCell.Offset(itime, 3) = temp(31, 2)

      ActiveCell.Offset(itime, 4) = temp(41, 2)

      ActiveCell.Offset(itime, 5) = temp(51, 2)

      ActiveCell.Offset(itime, 6) = temp1(1)

      ActiveCell.Offset(itime, 7) = temp1(2)

      ActiveCell.Offset(itime, 8) = temp1(3)

      ActiveCell.Offset(itime, 9) = temp1(4)

      ActiveCell.Offset(itime, 10) = temp1(5)

      Next itime

      Sheets("Sheet2").Select

      ActiveSheet.[A1].Select

      ActiveCell.Offset(0, 0) = "X"

      ActiveCell.Offset(0, 1) = "FDM:Fo=0.00125"

      ActiveCell.Offset(0, 2) = "FDM:Fo=0.0125"

      ActiveCell.Offset(0, 3) = "FDM:Fo=0.0625"

      ActiveCell.Offset(0, 4) = "FDM:Fo=0.125"

      ActiveCell.Offset(0, 5) = "FDM:Fo=0.25"

      ActiveCell.Offset(0, 6) = "theory:Fo=0.00125"

      ActiveCell.Offset(0, 7) = "theory:Fo=0.0125"

      ActiveCell.Offset(0, 8) = "theory:Fo=0.0625"

      ActiveCell.Offset(0, 9) = "theory:Fo=0.125"

      ActiveCell.Offset(0, 10) = "theory:Fo=0.25"

      For i = 1 To 101

      resultemp(i, 1) = dx * (i - 1)

      Next i

      Call theory_solution_space(resultemp, Fo, ramda, c, roh, SL)

      For i = 1 To 101

      ActiveCell.Offset(i, 0) = resultemp(i, 1)

      ActiveCell.Offset(i, 1) = resultemp(i, 2)

      ActiveCell.Offset(i, 2) = resultemp(i, 3)

      ActiveCell.Offset(i, 3) = resultemp(i, 4)

      ActiveCell.Offset(i, 4) = resultemp(i, 5)

      ActiveCell.Offset(i, 5) = resultemp(i, 6)

      ActiveCell.Offset(i, 6) = resultemp(i, 7)

      ActiveCell.Offset(i, 7) = resultemp(i, 8)

      ActiveCell.Offset(i, 8) = resultemp(i, 9)

      ActiveCell.Offset(i, 9) = resultemp(i, 10)

      ActiveCell.Offset(i, 10) = resultemp(i, 11)

      Next i

      End Sub

      Sub theory_solution_space(resultemp, Fo, ramda, hc, roh, SL)

      pai = 3.14159265358979

      AKAPPA = ramda / hc / roh / SL / SL

      For ii = 1 To 5

      t = Fo(ii) / AKAPPA

      For j = 1 To 101

      x = resultemp(j, 1)

      theta = 0#

      theta100 = 100#

      For i = 1 To 1000

      n = (i - 1) * 2 + 1

      an = n

      theta = theta + 1# / an * Exp(-an * an * pai * pai * ramda / hc / roh * t / SL / SL) * Sin(an * pai * x / SL)

      Next i

      resultemp(j, ii + 6) = theta100 * 4# / pai * theta

      Next j

      Next ii

      End Sub

      Sub theory_solution_time(time, temp1, ramda, c, roh, SL)

      pai = 3.14159265358979

      For j = 1 To 5

      x = 0.01 * j

      theta = 0#

      theta100 = 100#

      bkappa = ramda / c / roh / SL / SL

      For i = 1 To 1000

      n = (i - 1) * 2 + 1

      an = n

      theta = theta + 1# / an * Exp(-an * an * pai * pai * bkappa * time) * Sin(an * pai * x / SL)

      Next i

      temp1(j) = theta100 * 4# / pai * theta

      Next j

      End Sub

ここまでが、Excel　macro　FDM　です。

以下、FORTRANです。

FORTRAN 　SAMPLE 　PROGRAM

　（以下、Cの書いてある所から下をコピーしてTeraPadに

　　　貼り付け、 FDM_heat.forという名称で保存）

C

C      file name      FDM_heat.for

C

C   compile

C   >gfortran FDM_heat.for -o FDM_heat.exe

C

C

C   exequte

C   >FDM_heat

C

C

C-----------------------------------------------------------------

C   このプログラムを使うにはMinGWのパッケージが必要です。

C

C　　ネットで「MinGW」を検索してダウンロードする。

C　　インストールは、gfortran,C,C++のコンポーネントを入れる。

C

C　　compileと実行は上記の通り

C

C　　実行すると、FDM_space.txtとFDM_time_series.txtが出来る。

C

C　　FDM_space.txtは、Excelに貼り付け散布図を描くと、温度空間分布

C　　のグラフが描ける。

C

C　　FDM_time_series.txtもExcelに貼り付け散布図を描くと、時間軸の

C　　グラフが描ける。

C

C

C-----------------------------------------------------------------

C

      implicit real*8 (a-h,o-z)

      dimension temp(101,3),tempnew(101,3),resultemp(101,11)

      dimension Fo(5),IPTIME(5),temp1(5)

      dx=0.001

      dy=0.001

      dt=0.01

      SL=dx*100.0

      ramda=46.0

      c=460.5

      roh=7870.0

      Fo(1)=0.00125

      Fo(2)=0.0125

      Fo(3)=0.0625

      Fo(4)=0.125

      Fo(5)=0.25

      AKAPPA=ramda/c/roh

      do 10 i=1,5

      IPTIME(i)=Fo(i)/AKAPPA*SL*SL/dt

   10 continue

      A1=dt*ramda/c/roh/dx/dx

      A2=A1

      A3=dt*ramda/c/roh/dy/dy

      A4=A3

      AA=1.0-(A1+A2+A3+A4)

      OPEN(UNIT=8,FILE='FDM_time_series.txt',

     &STATUS='UNKNOWN')

      write(8,400)

  400 format("time,FDM:x=0.01m,FDM:x=0.02m,FDM:x=0.03m,FDM:x=0.04m,",

     &            "FDM:x=0.05m,",

     &            "theory:x=0.01m,theory:x=0.02m,theory:x=0.03m,",

     &            "theory:x=0.04m,theory:x=0.05m")

C

C

      do 35 i=1,101

      do 45 j=1,3

      temp(i,j)=100.0

   45 continue

   35 continue

      temp(1,1)=0.0

      temp(101,1)=0.0

      temp(1,2)=0.0

      temp(101,2)=0.0

      temp(1,3)=0.0

      temp(101,3)=0.0

C

C

      do 20 itime=1,IPTIME(5)+3

      time=itime*dt

      write(6,*)'TIME=',time 

      j=1

      do 30 i=2,100

      tempnew(i,j)=A1*temp(i+1,j)+A2*temp(i-1,j)

     &            +A3*temp(i,j+1)+A4*temp(i,j+1)+AA*temp(i,j)

   30 continue

      j=2

      do 40 i=2,100

      tempnew(i,j)=A1*temp(i+1,j)+A2*temp(i-1,j)

     &            +A3*temp(i,j-1)+A4*temp(i,j+1)+AA*temp(i,j)

   40 continue

      j=3

      do 50 i=2,100

      tempnew(i,j)=A1*temp(i+1,j)+A2*temp(i-1,j)

     &            +A3*temp(i,j-1)+A4*temp(i,j-1)+AA*temp(i,j)

   50 continue

      do 60 i=1,5

      if(IPTIME(i).eq.itime)then

      do 70 k=1,101

      resultemp(k,i+1)=tempnew(k,2)

   70 continue

      else

      endif

   60 continue

      do 80 i=1,101

      do 90 j=1,3

      temp(i,j)=tempnew(i,j)

   90 continue

   80 continue

      call theory_solution_time(time,temp1,5,ramda,c,roh,SL)

      write(8,300)time,temp(11,2),temp(21,2),temp(31,2),

     &                 temp(41,2),temp(51,2),(temp1(i),i=1,5)

   20 continue

C

C

C

      OPEN(UNIT=3,FILE='FDM_space.txt',

     &STATUS='UNKNOWN')

C

C

C

      write(3,200)

  200 format('X,FDM:Fo=0.00125,FDM:Fo=0.0125,FDM:Fo=0.0625,'

     &'FDM:Fo=0.125,FDM:Fo=0.25,'

     &'theory:Fo=0.00125,theory:Fo=0.0125,theory:Fo=0.0625,'

     &'theory:Fo=0.125,theory:Fo=0.25')

      do 15 i=1,101

      resultemp(i,1)=dx*float(i-1)

   15 continue

C

      call theory_solution_space(resultemp,101,11,Fo,5,

     &ramda,c,roh,SL)

C

      do 25 i=1,101

      write(3,300)(resultemp(i,j),j=1,11)

  300 format(E14.6,10(',',E14.6))

   25 continue

      CLOSE(3)

      stop

      end

C------------------------------------------------------------------

C

      subroutine theory_solution_space(resultemp,n1,n2,Fo,n3,

     &ramda,hc,roh,sl)

      implicit real*8 (a-h,o-z)

      dimension resultemp(n1,n2),Fo(n3)

      PAI=3.14159265358979

      AKAPPA=ramda/hc/roh/SL/SL

      do 30 ii=1,5

      t=Fo(ii)/AKAPPA

      do 20 j=1,101

      x=resultemp(j,1)

      theta=0.

      theta100=100.

      do 10 i=1,1000

      n=(i-1)*2+1

      an=float(n)

      theta=theta+1./an*EXP(-an*an*pai*pai*ramda/hc/roh*t/sl/sl)

     &      *sin(an*pai*x/sl)

   10 continue

      resultemp(j,ii+6)=theta100*4./pai*theta

   20 continue

   30 continue

      return

      end

C------------------------------------------------------------------

C

      subroutine theory_solution_time(time,temp1,n1,ramda,c,roh,SL)

      implicit real*8 (a-h,o-z)

      dimension temp1(n1)

      pai=3.14159265358979

      do 20 j=1,5

      x=0.01*float(j)

      theta=0.

      theta100=100.

      bkappa=ramda/c/roh/SL/SL

      do 10 i=1,1000

      n=(i-1)*2+1

      an=float(n)

      theta=theta+1./an*EXP(-an*an*pai*pai*bkappa*time)

     &      *sin(an*pai*x/SL)

   10 continue

      temp1(j)=theta100*4./pai*theta

   20 continue

      return

      end

計算結果　作図方法　Excel

実行するとFDM_space.txtができます。

開き、CTRL+A、CTRL+CですべてコピーしてExcelの[A1]セルに

貼り付けます。

すべての列を選択して散布図を描くとつぎのグラフが描けます。

Excelで時間軸のグラグは、FDM_time_series.txtを開き、

CTRL+A、CTRL+Cですべて選択して散布図を描けば得られます。

次に示します。縦軸（温度）をリニアと対数の両方を示します。