1
$\begingroup$

I have to write a solver for 2D equation: $$\partial_t u = u^2(\partial_x ^2 u + \partial_y ^2 u)$$ I try to use explicit method: $$\partial_t u = \frac{u_{i,j}^{k+1} - u_{i,j}^k}{\tau}$$ and $$\partial_x ^2 u = \frac{u_{i-1,j}^k - 2u_{i,j}^k + u_{i+1,j}^k}{h^2}$$ and tha same for $\partial_y ^2 u$ with points $j-1, j, j+1$. Also I use: $$ u_{i,j}^k = \frac{u_{i+1,j}^k + u_{i,j}^k}{2}$$ for the coefficient before $\partial_x ^2$ and the same for $y$ with points $j+1, j$. I wrote the following code in Mathematica:

nn = 6;
mm = 200;
h = 1./nn;
tau = 1./mm;
v0[x_, y_] := 1.;


Table[{x[i] = i*h, y[j] = j*h}, {i, 0, nn}, {j, 0, nn}];
Table[t[i] = i*t, {i, 0, mm}];
Table[u[i, j, 0] = v0[x[i], y[j]], {i, 0, nn}, {j, 0, nn}];
Table[u[0, j, k] = 0, {j, 0, nn - 1}, {k, 0, mm}];
Table[u[nn, j, k] = 0, {j, 0, nn - 1}, {k, 0, mm}];
Table[u[i, 0, k] = 0, {i, 0, nn - 1}, {k, 0, mm}];
Table[u[i, nn, k] = 0, {i, 0, nn - 1}, {k, 0, mm}];

Do[ 
  u[i, j, k + 1] = 
   u[i, j, k] + 
    tau*(((u[i - 1, j, k] - 2*u[i, j, k] + u[i + 1, j, k])/
          h^2)*(( u[i, j, k] + u[i - 1, j, k])/
           2)^2 + ((u[i, j - 1, k] - 2*u[i, j, k] + u[i, j + 1, k])/
          h^2)*((u[i, j, k] + u[i, j - 1, k])/2)^2), {k, 0, mm}, {i, 
   1, nn - 1}, {j, 1, nn - 1}];



g[k_] := Table[u[i, j, k], {i, 0, nn}, {j, 0, nn}];
p = Table[
  ListPlot3D[g[k], PlotRange -> {{0, nn}, {0, nn}, {0, 6}}], {k, 0, 
   mm, 20}]

I use the simpliest conditions: $u(0,0,t) = u(10,10,t) = 0$ and $u(x,y,0) = 1$

Improve this question
$\endgroup$
  • $\begingroup$ Maybe I'm missing something, but why is u interpolated along x and y for the coefficient of the u diffusion? I assume "t" in the second equation is supposed to be delta t. Also, please format the code you've provided and specify what language is being used. $\endgroup$ – Charles Apr 27 '17 at 2:49
  • 1
    $\begingroup$ What are the initial/boundary conditions? $\endgroup$ – zhk Apr 27 '17 at 8:59
  • $\begingroup$ I see a u^h^2 term in the Mathematica code that is almost certainly a bug. $\endgroup$ – Spencer Bryngelson Apr 27 '17 at 9:00
  • $\begingroup$ I've fixed the bug. Also I tried to use the coefficient as $u_{i,j,k}^2$ but it still doesn't work $\endgroup$ – Daniel Alexsandrovich Apr 27 '17 at 16:33
  • 2
    $\begingroup$ You might clarify what you mean by "it doesn't work". $\endgroup$ – Spencer Bryngelson Apr 27 '17 at 18:08