forked from burakbayramli/books
-
Notifications
You must be signed in to change notification settings - Fork 0
/
implicit.py
66 lines (55 loc) · 1.37 KB
/
implicit.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
import math, pylab
# fixed parameters
T0 = 25.0 # temperature gradient, C
D = 1e-4 # thermal diffusivity, m^2/s
L = 1.0 # length of rod, m
# input parameters
dx = L*input('grid spacing in units of rod length (L) -> ')
dt = (L**2/D)*input('time step in units of (L^2/D) -> ')
tmax = (L**2/D)*input('evolution time in units of (L^2/D) -> ')
# coefficients of the tridiagonal matrix
alpha = gamma = -D*dt/dx**2
beta = 1.0+2.0*D*dt/dx**2
# construct initial data
N = int(L/dx)
x = [0.0]*(N+1)
v = [0.0]*(N+1)
u = [0.0]*(N+1)
for j in range(N+1):
x[j] = j*dx
u[0] = T0
# prepare animated plot
pylab.ion()
(line, ) = pylab.plot(x, u, '-k')
pylab.ylim(0, T0)
pylab.xlabel('x (m)')
pylab.ylabel('Temperature (Celcius)')
# preform the evolution
t = 0.0
while t < tmax:
# update plot
line.set_ydata(u)
pylab.title('t = %5f'%t)
pylab.draw()
pylab.pause(0.1)
# swap u and v
(u, v) = (v, u)
# boundary conditions
u[0] = T0
u[N] = 0.0
# set the j=1 and j=N-1 points of v to the correct values
v[1] -= alpha*u[0]
v[N-1] -= gamma*u[N]
# forward sweep
u[1] = v[1]/beta
v[1] = gamma/beta
for j in range(2, N):
den = beta-alpha*v[j-1]
u[j] = (v[j]-alpha*u[j-1])/den
v[j] = gamma/den
# backward sweep
for j in reversed(range(1, N-1)):
u[j] -= u[j+1]*v[j]
t += dt
pylab.ioff()
pylab.show()