One dimensional Heat Equation Solving Algorithm

We let u(t, x) be the solution function for the differential equation defined by

• u_{t}(t, x) = α * u_{xx}(t, x) + S(t, x), 0 < x <= L, 0 < t <= T
• u(t, 0) = Φ_1(t) (Left Dirichlet Boundary Condition)
• u_{t}(t, L) = Φ_2(t) (Right Neumann Boundary Condition)
• u(0, x) = f(x) (initial values)

Let u[i][j] = u(i * dt, j * dx) for i = 0 to i = K - 1 and j = 0 to j = N - 1 for dt = T / (K - 1) and dx = L / (N - 1). For dt small enough, we can approximate u_{t}(t, x) with the forward difference:

u_{t}(t, x) = (u(t + dt, x) - u(t, x)) / dt

For dx small enough, we can approximate u_{xx}(t, x) with the second order central difference:

u_{xx}(t, x) = (u(t, x + dx) - 2 * u(t, x) + u(t, x - dx)) / (dx ** 2)

So after substituting and rearranging u_{xx}, u_{t} in u_{t}(t, x) = α * u_{xx}(t, x) + S(t, x), we get

u(t + dt, x) = u(t, x) + (α * dt / dx ** 2) * (u(t, x + dx) - 2 * u(t, x) + u(t, x - dx)) + S(t, x) * dt

Thus, we can compute all values of u[i][j] with the following:

u[i][j] = { u[i - 1, j] + (α * dt / dx ** 2) * (u[i - 1, j + 1] - 2 * u[i - 1, j] + u[i - 1, j - 1]) + S(i * dt, j * dx) * dt   0 < i <= K - 1, 0 < j < N - 1
          { Φ_1(i * dt)                                                                                                         j = N - 1
          { u[i, j - 2] + Φ_2(i * dt) * dt                                                                                      j = 0
          { f(j * dx)                                                                                                           i = 0

[full code implementation]