# How to discretize a nonlinear PDE fast?

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?

-
Would you format your code, please? – rcollyer Sep 11 '12 at 13:11
Your question is unreadable. Would you please consider to read the FAQ of this site first? Furthermore, please see here how you can format your question properly: mark code, use LaTeX, etc.. – halirutan Sep 11 '12 at 13:35
I revised the question, now please give some comments. – Fazlollah Soleymani Sep 11 '12 at 14:47
I looked at equations world web site, looked at all the non-linear parabolic pde's there, and do not see one that has this form you have? May be you can check this site and let us know which one of the PDE's listed in yours. Here is the link eqworld.ipmnet.ru/en/solutions/npde/npde-toc1.htm Either, when all else fails, you can start with finite difference method and go from there. – Nasser Sep 11 '12 at 17:21
it should be the nonlinear burgers equations Ok. But the Burger pde form that I know is not what you wrote. The Burgers viscous PDE is u_tt + u * u_x = mu * u_xx while what you wrote above is u_tt + u * u_xx = mu * u_xx i.e. u is multiplied by first derivative of u not the second partial derivative as you have there? Here is the wiki page en.wikipedia.org/wiki/Burgers%27_equation – Nasser Sep 11 '12 at 19:39

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,
"u(x,t) at t=" <>ToString[\[CapitalDelta]t n // N]]
],
{n, 0, timesteps,30}];


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.

-
 Whee, poles... :D – J. M.♦ Sep 12 '12 at 2:54 You solved Subscript[u, t] + u Subscript[u, xx] = c Subscript[u, xx]. Please consider solving the following NPDE instead: Subscript[u, t] + u Subscript[u, x] = c Subscript[u, xx]. – Fazlollah Soleymani Sep 12 '12 at 7:25 Sorry, I have to cast my first downvote…This solution is wrong. It just gives me the warning message \$RecursionLimit::reclim.I set some effort in correcting the mistake but sadly I failed, the only thing I know is that the definition of u[i,n] is incorrect.(See here for more details.) And, even with the solution in the link I give, the code won't work, I guess the reason is the recursion in this case is quite complex and it doesn't work as our expect. – xzczd Oct 3 '12 at 12:38