求助：一维非定常热传导方程BTCS(隐式格式)的求解Fortran程序

aliu

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

cora

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. 
   
8.       open(unit=9,file='btcs.out')
   
9. 
   
10. c...Read initial data samples.  Samples evenly spaced.
    
11. c...Data assumed periodic.
    
12.       open(unit=8,file='nb.dat',status='old')
    
13.       read(8,*) n, lambda, tfinal
    
14.       if(n.gt.nmax) then
    
15.         write(9,*) '****Too many data points****'
    
16.         close(unit=8)
    
17.         close(unit=9)
    
18.         stop
    
19.       endif
    
20.       if(n.lt.2) then
    
21.         write(9,*) '****Too few data points****'
    
22.         close(unit=8)
    
23.         close(unit=9)
    
24.         stop
    
25.       endif
    
26.       if(lambda.lt.0.01) then
    
27.         write(9,*) '****Lambda small or negative****'
    
28.         close(unit=8)
    
29.         close(unit=9)
    
30.         stop
    
31.       endif
    
32.       i=1
    
33.       read(8,*,err=1000,end=1000) xmin, u(1)
    
34.       do 10, i=2,n
    
35.         read(8,*,err=1000,end=1000) dummy, u(i)
    
36. 10   continue
    
37.       i=n+1
    
38.       read(8,*,err=1000,end=1000) xmax, u(n+1)
    
39.       if(abs(u(n+1)-u(1)).gt..0001) then
    
40.         write(9,*) '****Data not periodic****'
    
41.         close(unit=8)
    
42.         close(unit=9)
    
43.         stop
    
44.       endif
    
45.       if(xmax.le.xmin+.0001) then
    
46.         write(9,*) '****Bad x-axis****'
    
47.         close(unit=8)
    
48.         close(unit=9)
    
49.         stop
    
50.       endif
    
51.       u(0) = u(n)
    
52.       u(-1) = u(n-1)
    
53.       u(-2) = u(n-2)
    
54.       u(n+2) = u(2)
    
55.       do 15, i=1,n+1
    
56.         u0(i) = u(i)
    
57. 15   continue
    
58. 
    
59.       delta_x=(xmax-xmin)/real(n)
    
60.       delta_t=lambda*delta_x
    
61.       itert=nint(tfinal/delta_t)
    
62.       write(9,*) 'Final time requested: ', tfinal
    
63.       tfinal = real(itert)*delta_t
    
64.       write(9,*) 'Actual final time: ', tfinal
    
65.       write(9,*) 'delta_t = ', delta_t
    
66.       write(9,*) 'delta_x = ', delta_x
    
67.       write(9,*) 'lambda = ', lambda
    
68. 
    
69.       do 500, it=1,itert
    
70. 
    
71. c...Form ``periodic tridiagonal'' system of equations
    
72. 
    
73. c...BTCS
    
74.         do 20, i= 1,n
    
75. c...Lower sub-diagonal
    
76.           a(i) = -.5*lambda
    
77. c...Diagonal
    
78.           b(i) = 1.0
    
79. c...Upper sub-diagonal
    
80.           c(i) = .5 * lambda
    
81. c...Right hand side
    
82.           d(i) = u(i)
    
83. c...Upper right element
    
84.           e = -.5* lambda
    
85. c...Lower left element
    
86.           f = .5 * lambda
    
87. 20     continue
    
88. 
    
89. c...BTBS
    
90. c       do 20, i= 1,n
    
91. c...Lower sub-diagonal
    
92. c         a(i) = -lambda
    
93. c...Diagonal
    
94. c         b(i) = 1.0 + lambda
    
95. c...Upper sub-diagonal
    
96. c         c(i) = 0.
    
97. c...Right hand side
    
98. c         d(i) = u(i)
    
99. c...Upper right element
    
100. c         e = - lambda
     
101. c...Lower left element
     
102. c         f = 0.
     
103. c20     continue
     
104. 
     
105.         call ptrid(n,a,b,c,d,e,f,h)
     
106. 
     
107.         do 150, i=1,n
     
108.           u(i) = h(i)
     
109. 150    continue
     
110. 
     
111.         u(0) = u(n)
     
112.         u(-1) = u(n-1)
     
113.         u(-2) = u(n-2)
     
114.         u(n+1) = u(1)
     
115.         u(n+2) = u(2)
     
116. 
     
117. 500  continue
     
118. 
     
119.       write(9,*)
     
120.       write(9,1050)
     
121.       do 800, i=1,n+1
     
122.         write(9,1100) i, u0(i), u(i), abs(u(i)-u0(i))
     
123. 800  continue
     
124. 
     
125.       close(unit=8)
     
126.       close(unit=9)
     
127. 
     
128. C...write simple file for plotting
     
129.       open(unit=10,file='btcs.plt')
     
130.       do 900, i=1,n+1
     
131.         write(10,1150) -1 + real(2*i-2)/real(n), u(i),u0(i)
     
132. 900  continue
     
133.       close(unit=10)
     
134. 
     
135.       stop
     
136. 
     
137. 1000 write(9,*) '****Error reading data point number ', i,'****'
     
138.       close(unit=8)
     
139.       close(unit=9)
     
140.       stop
     
141. 
     
142. 1050 format(4x,'N',11x,'INITIAL',12x,'BTCS',11x, 'DIFFERENCE')
     
143. 1100 format(i5,5x,f13.8,5x,f13.8,5x,f13.8)
     
144. 1150 format(f15.9,5x,e25.9,5x,f15.9)
     
145. 
     
146.       end
     
147. 
     
148.       subroutine ptrid(n,a,b,c,d,e,f,x)
     
149.       integer n
     
150.       double precision a(n),b(n),c(n),d(n),e,f,x(n)
     
151.       double precision r(n),s(n)
     
152. 
     
153.       r(1) = e
     
154.       s(1) = f
     
155. 
     
156. c...Zero out the subdiagonal
     
157.       do 10, i=2,n-2
     
158. c       a(i) = a(i) - a(i)*b(i-1)/b(i-1)
     
159.         b(i) = b(i) - a(i)*c(i-1)/b(i-1)
     
160.         d(i) = d(i) - a(i)*d(i-1)/b(i-1)
     
161.         r(i) =      - a(i)*r(i-1)/b(i-1)
     
162. 10   continue
     
163. 
     
164. c     a(n-1) = a(n-1) - a(n-1)*b(n-2)/b(n-2)
     
165.       b(n-1) = b(n-1) - a(n-1)*c(n-2)/b(n-2)
     
166.       c(n-1) = c(n-1) - a(n-1)*r(n-2)/b(n-2)
     
167.       d(n-1) = d(n-1) - a(n-1)*d(n-2)/b(n-2)
     
168. 
     
169. c     a(n) = a(n) - a(n)*b(n-1)/b(n-1)
     
170.       b(n) = b(n) - a(n)*c(n-1)/b(n-1)
     
171.       d(n) = d(n) - a(n)*d(n-1)/b(n-1)
     
172. 
     
173.       if(abs(f).gt..0001) then
     
174. c...Chase the non-zero element from left to right across the
     
175. c...bottom row.
     
176.         do 15, i=1,n-2
     
177. c         s(i)   =  s(i) - s(i)*b(i)/b(i)
     
178.           s(i+1) =       - s(i)*c(i)/b(i)
     
179.           b(n)   =  b(n) - s(i)*r(i)/b(i)
     
180.           d(n)   =  d(n) - s(i)*d(i)/b(i)
     
181. 15     continue
     
182. c         s(n-1) =  s(n-1) - s(n-1)*b(n-1)/b(n-1)
     
183.           b(n)   =  b(n)   - s(n-1)*c(n-1)/b(n-1)
     
184.           d(n)   =  d(n)   - s(n-1)*d(n-1)/b(n-1)
     
185.       endif
     
186. 
     
187. c...Solve the upper-triangular system using
     
188. c...back-substitution.
     
189.       x(n) = d(n)/b(n)
     
190.       x(n-1) = (d(n-1)-c(n-1)*x(n))/b(n-1)
     
191.       do 20, i=n-2,1,-1
     
192.         x(i) = (d(i)-c(i)*x(i+1)-r(i)*x(n))/b(i)
     
193. 20   continue
     
194. 
     
195.       return
     
196.       end

复制代码

aliu

感谢您的帮助！:handshake