声振论坛

 找回密码
 我要加入

QQ登录

只需一步,快速开始

查看: 2798|回复: 2

[基础理论] 求助:一维非定常热传导方程BTCS(隐式格式)的求解Fortran程序

[复制链接]
发表于 2007-11-19 22:50 | 显示全部楼层 |阅读模式

马上注册,结交更多好友,享用更多功能,让你轻松玩转社区。

您需要 登录 才可以下载或查看,没有账号?我要加入

x
我是个程序新手,但是现在急需要求解一维非定常热传导方程的隐式格式,需要永fortran或者
matlab的程序,盼望好心的高手帮帮忙啊!不生感激!
回复
分享到:

使用道具 举报

发表于 2007-11-21 09:43 | 显示全部楼层
  1.       program btcs
  2. c...Performs backward-time central-space  (IMPLICIT EULER)
  3.       parameter (nmax=100,delta=.000001)
  4.       double precision lambda,u(-2:nmax+2),h(1:nmax),u0(1:nmax+1)
  5.       double precision a(1:nmax), b(1:nmax), c(1:nmax)
  6.       double precision d(1:nmax), e, f

  7.       open(unit=9,file='btcs.out')

  8. c...Read initial data samples.  Samples evenly spaced.
  9. c...Data assumed periodic.
  10.       open(unit=8,file='nb.dat',status='old')
  11.       read(8,*) n, lambda, tfinal
  12.       if(n.gt.nmax) then
  13.         write(9,*) '****Too many data points****'
  14.         close(unit=8)
  15.         close(unit=9)
  16.         stop
  17.       endif
  18.       if(n.lt.2) then
  19.         write(9,*) '****Too few data points****'
  20.         close(unit=8)
  21.         close(unit=9)
  22.         stop
  23.       endif
  24.       if(lambda.lt.0.01) then
  25.         write(9,*) '****Lambda small or negative****'
  26.         close(unit=8)
  27.         close(unit=9)
  28.         stop
  29.       endif
  30.       i=1
  31.       read(8,*,err=1000,end=1000) xmin, u(1)
  32.       do 10, i=2,n
  33.         read(8,*,err=1000,end=1000) dummy, u(i)
  34. 10   continue
  35.       i=n+1
  36.       read(8,*,err=1000,end=1000) xmax, u(n+1)
  37.       if(abs(u(n+1)-u(1)).gt..0001) then
  38.         write(9,*) '****Data not periodic****'
  39.         close(unit=8)
  40.         close(unit=9)
  41.         stop
  42.       endif
  43.       if(xmax.le.xmin+.0001) then
  44.         write(9,*) '****Bad x-axis****'
  45.         close(unit=8)
  46.         close(unit=9)
  47.         stop
  48.       endif
  49.       u(0) = u(n)
  50.       u(-1) = u(n-1)
  51.       u(-2) = u(n-2)
  52.       u(n+2) = u(2)
  53.       do 15, i=1,n+1
  54.         u0(i) = u(i)
  55. 15   continue

  56.       delta_x=(xmax-xmin)/real(n)
  57.       delta_t=lambda*delta_x
  58.       itert=nint(tfinal/delta_t)
  59.       write(9,*) 'Final time requested: ', tfinal
  60.       tfinal = real(itert)*delta_t
  61.       write(9,*) 'Actual final time: ', tfinal
  62.       write(9,*) 'delta_t = ', delta_t
  63.       write(9,*) 'delta_x = ', delta_x
  64.       write(9,*) 'lambda = ', lambda

  65.       do 500, it=1,itert

  66. c...Form ``periodic tridiagonal'' system of equations

  67. c...BTCS
  68.         do 20, i= 1,n
  69. c...Lower sub-diagonal
  70.           a(i) = -.5*lambda
  71. c...Diagonal
  72.           b(i) = 1.0
  73. c...Upper sub-diagonal
  74.           c(i) = .5 * lambda
  75. c...Right hand side
  76.           d(i) = u(i)
  77. c...Upper right element
  78.           e = -.5* lambda
  79. c...Lower left element
  80.           f = .5 * lambda
  81. 20     continue

  82. c...BTBS
  83. c       do 20, i= 1,n
  84. c...Lower sub-diagonal
  85. c         a(i) = -lambda
  86. c...Diagonal
  87. c         b(i) = 1.0 + lambda
  88. c...Upper sub-diagonal
  89. c         c(i) = 0.
  90. c...Right hand side
  91. c         d(i) = u(i)
  92. c...Upper right element
  93. c         e = - lambda
  94. c...Lower left element
  95. c         f = 0.
  96. c20     continue

  97.         call ptrid(n,a,b,c,d,e,f,h)

  98.         do 150, i=1,n
  99.           u(i) = h(i)
  100. 150    continue

  101.         u(0) = u(n)
  102.         u(-1) = u(n-1)
  103.         u(-2) = u(n-2)
  104.         u(n+1) = u(1)
  105.         u(n+2) = u(2)

  106. 500  continue

  107.       write(9,*)
  108.       write(9,1050)
  109.       do 800, i=1,n+1
  110.         write(9,1100) i, u0(i), u(i), abs(u(i)-u0(i))
  111. 800  continue

  112.       close(unit=8)
  113.       close(unit=9)

  114. C...write simple file for plotting
  115.       open(unit=10,file='btcs.plt')
  116.       do 900, i=1,n+1
  117.         write(10,1150) -1 + real(2*i-2)/real(n), u(i),u0(i)
  118. 900  continue
  119.       close(unit=10)

  120.       stop

  121. 1000 write(9,*) '****Error reading data point number ', i,'****'
  122.       close(unit=8)
  123.       close(unit=9)
  124.       stop

  125. 1050 format(4x,'N',11x,'INITIAL',12x,'BTCS',11x, 'DIFFERENCE')
  126. 1100 format(i5,5x,f13.8,5x,f13.8,5x,f13.8)
  127. 1150 format(f15.9,5x,e25.9,5x,f15.9)

  128.       end

  129.       subroutine ptrid(n,a,b,c,d,e,f,x)
  130.       integer n
  131.       double precision a(n),b(n),c(n),d(n),e,f,x(n)
  132.       double precision r(n),s(n)

  133.       r(1) = e
  134.       s(1) = f

  135. c...Zero out the subdiagonal
  136.       do 10, i=2,n-2
  137. c       a(i) = a(i) - a(i)*b(i-1)/b(i-1)
  138.         b(i) = b(i) - a(i)*c(i-1)/b(i-1)
  139.         d(i) = d(i) - a(i)*d(i-1)/b(i-1)
  140.         r(i) =      - a(i)*r(i-1)/b(i-1)
  141. 10   continue

  142. c     a(n-1) = a(n-1) - a(n-1)*b(n-2)/b(n-2)
  143.       b(n-1) = b(n-1) - a(n-1)*c(n-2)/b(n-2)
  144.       c(n-1) = c(n-1) - a(n-1)*r(n-2)/b(n-2)
  145.       d(n-1) = d(n-1) - a(n-1)*d(n-2)/b(n-2)

  146. c     a(n) = a(n) - a(n)*b(n-1)/b(n-1)
  147.       b(n) = b(n) - a(n)*c(n-1)/b(n-1)
  148.       d(n) = d(n) - a(n)*d(n-1)/b(n-1)

  149.       if(abs(f).gt..0001) then
  150. c...Chase the non-zero element from left to right across the
  151. c...bottom row.
  152.         do 15, i=1,n-2
  153. c         s(i)   =  s(i) - s(i)*b(i)/b(i)
  154.           s(i+1) =       - s(i)*c(i)/b(i)
  155.           b(n)   =  b(n) - s(i)*r(i)/b(i)
  156.           d(n)   =  d(n) - s(i)*d(i)/b(i)
  157. 15     continue
  158. c         s(n-1) =  s(n-1) - s(n-1)*b(n-1)/b(n-1)
  159.           b(n)   =  b(n)   - s(n-1)*c(n-1)/b(n-1)
  160.           d(n)   =  d(n)   - s(n-1)*d(n-1)/b(n-1)
  161.       endif

  162. c...Solve the upper-triangular system using
  163. c...back-substitution.
  164.       x(n) = d(n)/b(n)
  165.       x(n-1) = (d(n-1)-c(n-1)*x(n))/b(n-1)
  166.       do 20, i=n-2,1,-1
  167.         x(i) = (d(i)-c(i)*x(i+1)-r(i)*x(n))/b(i)
  168. 20   continue

  169.       return
  170.       end
复制代码
 楼主| 发表于 2007-11-27 20:30 | 显示全部楼层
感谢您的帮助!:handshake
您需要登录后才可以回帖 登录 | 我要加入

本版积分规则

QQ|小黑屋|Archiver|手机版|联系我们|声振论坛

GMT+8, 2024-11-11 04:48 , Processed in 0.070161 second(s), 18 queries , Gzip On.

Powered by Discuz! X3.4

Copyright © 2001-2021, Tencent Cloud.

快速回复 返回顶部 返回列表