I am having problems with discretizing the solution to the 1-D Wave equation. The problem lies in defining the size of the matrix needed to hold the solution for $u(x,t)$ for $0\leq x \lt 4 \pi$ and $0 \leq t \leq 16 \pi$. I'll outline my work thus far: It's a little long winded to get to the bit where I'm having problems, so I'll try to be concise:
$$\partial_t^2u(x,t) = c^2\partial_x^2u(x,t) $$
Using substitution to get two first order in time PDE's:
$$ \\ \partial_tu(x,t) = v(x,t) \\ \partial_t v(x,t) = c^2\partial_x^2u(x,t) $$
Adding a parabolic term to give stability to the numeric method gives:
$$ \\ \partial_t u(x,t) = v(x,t) + a\tilde{\partial_x^2}u(x,t) \\ \partial_t v(x,t) = c^2\partial_x^2u(x,t) + b\tilde{\partial_x^2}v(x,t) \\ $$
Applying the following finite difference operators to the above coupled first order equations,
$$ [\partial_t f(x,t)]_j^n = \frac{f_j^{n+1} - f^n_j}{\Delta t} \\ [\partial_x^2 f(x,t)]_j^n = \frac{f_{j+1}^{n} - 2f^n_j + f^n_{j-1}}{(\Delta x)^2} \\ [\tilde{\partial_x^2} f(x,t)]_j^n = \frac{f_{j+1}^{n} - 2f^n_j + f^n_{j-1}}{2\Delta t} $$
gives our numerical methods to solve u and v:
$$ u_j^{n+1} = u_j^n + \Delta t v_j^n + \frac{a}{2}(u_{j+1}^n - 2u_j^n + u_{j-1}^n) \\v_j^{n+1} = v_j^n + \frac{\Delta t c^2}{(\Delta x)^2}(u_{j-1}^n - 2u_j^n + u_{j-1}^n) + \frac{b}{2}(v_{j-1}^n - 2v_j^n + v_{j-1}^n) $$
Initial Conditions are given as: $$u(x,0) = cox(x)$$ and $$\partial_tu(x,0) = v = 0$$
This implies:
$$ u_j^0 = cos(x_0 + j\Delta x), \forall j$$
Finally, boundary conditions are given as:
$$u(x,t) = u(x + 4\pi, t)$$ and $$v(x,t) = \partial_t(x,t) = \partial_tu(x + 4\pi, t)$$
And now here is where I am banging my head against the wall. I need to figure out the size of the matrix I need for a given time/space grid. Seems easy enoug, right? Just use $\frac{4\pi}{\Delta x}$ columns (plus two for the boundary conditions) and do the same for the time! However, this leads to non integer values, and I need to index my matrix using integers. Further, when I'm all done I have to fish out of my matrix $u(\pi, t)$ and compare it with the actual solution of my PDE. Figuring out which j index holds the value closest to $x=\pi$ using this method is also frustrating me.
What I've done is use the ceiling of $\frac{4\pi}{\Delta x}$ and $\frac{16\pi}{\Delta t}$ to "round" up to the number of rows and columns needed, but this feels very, very bad. I'm looking for an alternative, but can't think of one.
I have coded all of the following in python with the following method (if it is pertinent):
Dt = (0.25*DeltaX)/c
# Dt = (16*pi)/3
rows = ceil((16*pi) / Dt)
cols = ceil((4*pi) / DeltaX + 3)
u = zeros((rows, cols))
v = zeros((rows, cols))
for j in range(cols):
u[0][j] = cos(-DeltaX + j*DeltaX)
v[0][j] = 0
for i in range(1, rows):
for j in range(1, cols-1):
# u(x,t): fill row (time) i of grid space for columns (space) j
u[i][j] = u[i-1][j] + Dt*(v[i-1][j]) + (a/2)*(u[i-1][j+1] - 2*u[i-1][j] + u[i-1][j-1])
# v(x,t): fill row (time) i of grid space for columns (space) j
v[i][j] = v[i-1][j] + ((Dt*pow(c, 2))/pow(DeltaX, 2))*(u[i-1][j+1] - 2*u[i-1][j] + u[i-1][j-1])
+ (b/2)*(v[i-1][j+1] - 2*v[i-1][j] + v[i-1][j-1])
u = apply_bcs(u, DeltaX, i)
v = apply_bcs(v, DeltaX, i)
print(u)