Finite Difference Method in Mathematica
Published:
I wrote this notes for introducing the finite difference method. Basic heat equation and 2d wave equation are selected as cases, and the related Mathematica program are shown.
Heat equation in a slab:
\[\frac{\partial}{\partial t}T(x,t)=\chi \frac{\partial^2}{\partial x^2}T(x,t)+S(x,t) \\ T(0,t)=T_1 (t), \quad T(L,t)=T_2 (t)\]FTCS (Forward Euler in Time and Central difference in Space)
\[\text{Let} \quad T(x_j,t_n)=T_j^n \\ \frac{\partial}{\partial t}T(x_j,t_n)=\frac{T_j^{n+1}-T_j^n}{\Delta t}\\ \frac{\partial^2}{\partial x^2}T(x_j,t_n)=\frac{T_{j+1}^n-2T_j^n+T_{j-1}^n}{\Delta x^2}\\ T_j^{n+1}=T_j^n+\Delta t S_j^n+\frac{\chi \Delta t}{\Delta x^2}(T_{j+1}^n-2T_j^n+T_{j-1}^n)\] T[j_, n_] :=
T[j, n] =
T[j, n - 1] + \[CapitalDelta]t S[j,
n - 1] + \[Chi] \[CapitalDelta]t/\[CapitalDelta]x^2 (T[j + 1,
n - 1]
- 2 T[j, n - 1] + T[j - 1, n - 1])
T[0, n_] := T1[t] /. t -> n \[CapitalDelta]t
M = 10; T[M, n_] := T2[t] /. t -> n \[CapitalDelta]t
T[j_, 0] := T0[x] /. x -> j \[CapitalDelta]x
L = 1; \[Chi] = 1/4; \[CapitalDelta]t = 0.01; \[CapitalDelta]x = L/M;
S[j_, n_] := 0;
T0[x_] := x^2 (1 - x);
T1[t_] := t/(1 + 3 t);
T2[t_] := 0;
Table[ListPlot[Table[{\[CapitalDelta]x j, T[j, n]}, {j, 0, M}],
Joined -> True,
PlotRange -> {0, .4}, AxesLabel -> Automatic], {n, 0, 80, 20}]
2D wave equation
Mx = 30; My = 30;
Lx = Ly = 1;
\[CapitalDelta]x = Lx/Mx; \[CapitalDelta]y = Ly/My;
\[CapitalDelta]t = 0.2;
x[j_] := j \[CapitalDelta]x;
y[k_] := k \[CapitalDelta]y;
z0[j_, k_] := 0;
v0[j_, k_] := 0;
z[j_, k_, 0] := z0[j, k];
cw[x_, y_] := 1/10 (1 - .7 Exp[-2 (x - 0.5)^2 - 5 (y - 0.5)^2]);
c[j_, k_] := cw[x[j], y[k]];
z[0, k_, n_] := 0;
z[Mx, k_, n_] := 0;
z[j_, 0, n_] := x[j] (1 - x[j]) Sin[2 n \[CapitalDelta]t];
z[j_, My, n_] := 0;
z[j_, k_, n_] :=
(z[j, k, n] = 2 z[j, k, n - 1] - z[j, k, n - 2] +
\[CapitalDelta]t^2 c[j, k]^2/\[CapitalDelta]x^2
((z[j + 1, k, n - 1] - 2 z[j, k, n - 1] +
z[j - 1, k, n - 1]) + (z[j, k + 1, n - 1] -
2 z[j, k, n - 1] + z[k, k - 1, n - 1]))) /; n >= 2
z[j_, k_, 1] := 0;
zsol[n_] :=
zsol[n] =
Interpolation[
Flatten[Table[{j \[CapitalDelta]x, k \[CapitalDelta]y,
z[j, k, n]}, {j, 0, Mx}, {k, 0, My}], 1]]
Table[Plot3D[Evaluate[zsol[n][x, y]], {x, 0, 1}, {y, 0, 1},
PlotPoints -> 60, PlotRange -> Full, Mesh -> None], {n, 0, 140, 25}]
