How to discretize a nonlinear PDE fast?

Question

I wish to numerically solve the following PDE. Although there are some complete discussions for solving PDEs in tutorial/NDSolvePDE, there is no hint for the nonlinear case by discretization. Thus, I will be thankful to receive some helps on the following NPDE where $x \in [0,1]$, $t \in [0,2]$,

M=8; NN=8; m = M - 1; n = NN - 1; alpha = 5.; 
beta = 4.; c = 0.05; T = 2.; h = (1. - 0.)/M; 
k =T/NN

The problem is Subscript[u, t] + u Subscript[u, x] = c Subscript[u, xx] with initial condition

Subscript[u, x, 0] = (2. c beta Pi Sin[Pi Subscript[x, i]])
                   / (alpha + beta Cos[Pi Subscript[x, i]]) 

and boundary condition

Subscript[u, 0, t] = 0; Subscript[u, 1, t] = 0;

I tried the backward finite difference (FD) for Subscript[u, t] and the central FD for the others. I wrote the following code, but I think there are some gaps in it. Because, the approximate solutions do not match the exact one

u[x_, t_] := (2 c beta Pi Exp[-c Pi^2 t] Sin[Pi x])
            / (alpha + beta Exp[-c Pi^2 t] Cos[Pi x]);

Note that Subscript[w,i,j] stands for the approximation in the grid point $(x_i,t_j)$.

M = 8; NN = 8; m = M - 1; 
n = NN - 1; alpha = 5.; 
beta = 4.; c = 0.05; 
T = 2.; h = (1. - 0.)/M; 
k =T/NN;

(*Defining the Grid points*)
Table[Subscript[x, i] = 0 + i h, {i, 0, M}]; 
Table[Subscript[t, j] = 0 + j k, {j, 0, NN}];

(*Defining the Initial Conditions*)
For[i = 1, i <= m, i++, 
  Subscript[w, i, 0] = (2. c beta Pi Sin[Pi Subscript[x, i]])
                     / (alpha + beta Cos[Pi Subscript[x, i]])
];

(*Defining the Boundary Conditions*)
For[j = 1, j <= n, j++, 
  Subscript[w, 0, j] = 0
]; 
For[j = 1, j <= n, j++, 
  Subscript[w, 1, j] = 0
];

(*Defining the nonlinear equations due to discretization*)
For[i = 1, i <= m, i++,
  {
   For[j = 1, j <= n, j++,
     f[i, j] = Subscript[w, i, j] 
       + (k/(2 h)) Subscript[w, i, j] (Subscript[w, i + 1, j] - Subscript[w, i - 1, j]) 
       - (c k/(h^2)) (Subscript[w, i + 1, j] - 2 Subscript[w, i, j] 
                      + Subscript[w, i - 1, j]) 
       - Subscript[w,i, j - 1]
   ]
  }
];

F = Flatten[Table[f[i, j], {i, 1, m}, {j, 1, n}]]; 
Dimensions[F];
F // MatrixForm;
Vec = Flatten[Table[Subscript[w, i, j], {i, 2, M}, {j, 1, n}]];

(*Finding the solutions*)
Sol = Part[NSolve[F, Vec, Reals], 1]

Any suggestion is appreciated. In fact, what would be the final nonlinear system of equations resulting of discretization?

Answer 1

A quick answer now. I will come back to this once I have more time. First we use the common finite difference operators to discretize PDE.

Symbolics:

Clear[u];
RFSDiscrit[eq_] := 
Module[{mid},
       mid = Distribute@FullSimplify@ExpandAll[(
        eq /.{
          u[x, t] -> u[i, n],
          D[u[x, t], t] -> (u[i, n] - u[i, n - 1])/\[CapitalDelta]t,
          D[u[x, t], x, x] ->
          (u[i + 1, n - 1] - 2 u[i, n - 1] + u[i - 1, n - 1])/\[CapitalDelta]x^2
             }
        )];
       (First@Solve[mid == 0, u[i, n]])[[All, 2]] // First
       ];

Recursion:

Lets generate the recursive formula for the main equation.

exp = RFSDiscrit[D[u[x, t], t] + (c - u[x, t]) D[u[x, t], x, x]]

Now comes the rest of the recursion definition part.

alpha = 5.;
beta = 4.;
c = .05;
(*Use recursive formula*)
u[i_, n_] := u[i, n] = exp;
u[0, n_] := u1[t] /. t -> n \[CapitalDelta]t;
(*Number of discretization points for [0,L]*)
spacediscrit = 100; 
u[spacediscrit, n_] := u2[t] /. t -> n \[CapitalDelta]t;
u[i_, 0] := u0[x] /. x -> i \[CapitalDelta]x;
(* End of spatial boundary*)
L = 1;
timesteps = 1000;
\[CapitalDelta]t = 2/timesteps;
\[CapitalDelta]x = L/spacediscrit;
(*Initial Condition*)
u0[x_] = (2. c beta Pi Sin[Pi x])/(alpha + beta Cos[Pi x]);
(*Boundary Condition*)
u1[t_] = 0;
u2[t_] = 0;
(* Time snapshots for u(x,t_n)*)
pics = Table[ListPlot[Table[{\[CapitalDelta]x i, u[i, n]},{i, 0, spacediscrit}], 
            InterpolationOrder -> 3, Joined -> True, Mesh -> 10, 
            MeshStyle -> Red, Filling -> Axis, FillingStyle -> LightPink, 
            Frame -> True, PlotRange -> All, Axes -> None, 
            PlotLabel ->TraditionalForm[
             "u(x,t) at t=" <>ToString[\[CapitalDelta]t n // N]]
             ], 
            {n, 0, timesteps,30}];

Good/Bad-ness of...!

Now we can get $u(x,t)$ as a continuous function as follows.

dat = Table[Table[N@{\[CapitalDelta]x j, \[CapitalDelta]t n, u[j, n]}, 
            {j, 0,spacediscrit}],
      {n, 0, timesteps}];
fun = Interpolation[Flatten[dat, 1]]

InterpolatingFunction[{{0.,1.},{0.,2.}},<>]

Now a nice plot to celebrate the nasty numerical instability of the above equation on the given parameter setting.

The MMA function NDSolve FiniteDifferenceDerivative can also be used for discretization.