Numerical solution of nonlinear thermoconductivity equation

Question

I have to find and plot a numerical solution for the following equation (I have to write a solver):

$$u_{t} = (u^2 u_x)_x$$ with the following conditions $u(0,t) =0, u(1,t) = \sqrt{\frac{2c-2}{t}}, u(x,0) = 0$, where $c$ is a speed of the wave

I try to use the following method:

$$\frac{u_{m}^{n+1} - u_{m}^n}{\tau} = \frac{F_{m +1/2}^n - F_{m-1/2}^n}{h}$$ where $$F_{m+1/2}^n = (u_{m+1/2}^n)^2(\frac{u_{m+1}^n - u_{m-1}^n}{h})$$ where $$u_{m+1/2}^n = \frac{u_{m+1}^n + u_{m}^n}{2}$$

tt = 0.1;

tnn = 10000;

xnn = 20;

tau = tt/tnn;

h = 1/xnn;

u0[x_] := x;

Do[v[m, 0] = u0[h*m], {m, 0, xnn}]

Do[v[0, n] = 0; v[xnn, n] = 0;
 Do[v[m, n + 1] = 

v[m, n] + tau/h*( ((v[m + 1, n] - v[m, n])/ h)(((v[m + 1, n])^2 + (v[m, n])^2)/ 2) - ((v[m, n] - v[m - 1, n])/ h)(((v[m, n])^2 + (v[m - 1, n])^2)/2) ), {m, 1, xnn - 1}], {n, 0, tnn - 1}]

u = ListInterpolation[
   Table[v[i, j], {i, 0, xnn}, {j, 0, tnn}], {{0, 1}, {0, tt}}];

Plot3D[u[x, t], {x, 0, 1}, {t, 0, tt}, PlotRange -> All]

I have changed the code and have got the correct solvation, but for the u(0,t)=u(1,t) = 0, u(x,0) = f(x). I don't actually understand how should I write the conditions like u(1,t) = f(t).

Answer 1

You got good advices in comments, but as it is only 1D problem, you can manage it using fully explicit method with very small time steps. I checked your code in Mathematica and there are several misunderstandings there. After correcting them, the code works for me.

Similarly, your description in the question contains, I suppose, some mistakes. As I guess that this is some student project, so let me only indicate what shall be corrected ;-).

Your method in the question can be viewed as finite volume method for computational cells $(x_{m-1/2},x_{m+1/2})$. The quantity $F_{m +1/2}^n$ approximates the flux at the right point, the approximation in your question is inappropriate, check it, it shall use so-called central difference at point $x_{m+1/2}$. Analogously the nonlinear conductivity coefficient $u_{m+1/2}^n$ at $x_{m+1/2}$ is approximated inproperly.

Your code anyway does not implement what you describe in your question. So, check the scheme and implement it exactly and I suppose your code can work.

Be careful about a correct input data for your model. In your implementation you use different boundary and initial conditions as you describe in the question.

I hope after these indications you manage to finish the task, good luck!