-
Notifications
You must be signed in to change notification settings - Fork 3
/
Copy pathstep3.py
109 lines (81 loc) · 2.44 KB
/
step3.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
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
# -*- coding: utf-8 -*-
"""
Created on Thu Jul 7 16:05:12 2016
@author: clay_budin
Step 3: 1D diffusion
du/dt - v*d^2u/dx^2 = 0
v is diffusion parameter (viscocity)
if const then exact solns exist of form: u = u0 * e^i(kx-wt) w = 2*pi*f - plane wave
u = u0 * e^(ikx-vk^2t) - exponentially damped wave (when v > 0)
if v < 0 have exponential growth -> explosion
Also called heat equation or Fourier eq
note sign is different from convection
Finite Difference approximations to d^2u/dx^2
Central diff: d^2u/dx^2 ~= (u(x+1) - 2*u(x) + u(x-1)) / (delta_x^2) O(delta_x^2)
Central diff in time (isotropic)
Backward diff in space
Domain: [0,2]
Range: [1,2]
Initial Conditions: square wave, half-sine, full inverted cosine
Boundary Conditions: u = 1 @ x = 0,2
Now v*dt/(dx*dx) appears to be a critical factor
< 1 to be stable
<< 1 to be smooth, eg .5 produces stairstep jitter, .25 better
"""
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.animation as animation
# graphing vars
fig, ax = plt.subplots()
line, = ax.plot([], [], lw=1)
ax.grid()
# simulation constants
nx = 201
nt = 150
dt = 0.01
dx = 2.0 / (nx-1.0)
v = .005
# simlation grids - 1D
xdata, ydata, yprev = [], [], []
# global var - current time step
ct = 0
def init ():
global ct
if ct == 0: print "init() dx = " + str(dx) + " dt = " + str(dt) + " v*dt/(dx*dx) = " + str(v*dt/(dx*dx))
ct = 0
ax.set_xlim(-.1, 2.1)
ax.set_ylim(0.0,2.5)
del xdata[:]
del ydata[:]
for i in range(nx):
x = i*dx
xdata.append(x)
y = 1.0
# initialize range [.5,1] with various test fns
if x >= .5 and x <= 1.0: y = 2.0 # square wave
#if x >= .5 and x <= 1.0: y = 1.0 + np.sin((x-.5)*np.pi*2.0) # half-sine
# if x >= .5 and x <= 1.0: # full inverted cosine
# #print "t = " + str(np.cos((x-.5)*np.pi*4.0))
# y = 1.0 + .5*(1.0 - np.cos((x-.5)*np.pi*4.0))
# #print "y = " + str(y)
#y -= .5
#if x >= .5 and x <= 1.0: y = 1.0 + .5*np.sin((x-.5)*np.pi*4.0) # full-sine
ydata.append(y)
yprev.append(y)
line.set_data(xdata, ydata)
return line,
def data_gen ():
global ct
while ct <= nt:
#print "data_gen() ct =", ct
ct += 1
for i in range(nx): yprev[i] = ydata[i]
for i in range(1,nx-1):
ydata[i] = yprev[i] + (v*dt/(dx*dx))*(yprev[i+1]-2.0*yprev[i]+yprev[i-1])
yield
def run (data):
#print "run()"
line.set_data(xdata, ydata)
return line,
ani = animation.FuncAnimation(fig, run, data_gen, blit=False, interval=10, repeat=True, init_func=init)
plt.show()