Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

I want to solve the following heat conduction equation using numerical methods:

D[u[x, t], t] -alpha*D[u[x, t], {x, 2}] == 0

u[x, 0] == 1/(1 + x^2)^0.25,

u[-10, t] == u[10, t] == 0,

x0 = -10;
xn = 10;
dx = 0.1;
n = 200;
t0 = 0;
tn = 10;
dt = 0.01;
steps = tn/dt;
alpha = 0.2;

I define:

points = Table[-10 + i*dx, {i, 0, n}];

u[0] = Module[{i, l= {0}},  For[i = 1, i <= n - 1, i = i + 1; 
AppendTo[l, 1/(1 + points[[i]]^2)^0.25]]; AppendTo[l, 0];l];    

previous[dt_, dx_, u0_, alpha_] := Module[{i, ans}, 
ans = u0[[2]] + alpha*(dt/dx^2)*(u0[[3]] - 2 u0[[2]] + u0[[1]]); ans];

I aproximate the values of the following steps:

u[k_] := Module[{l,ans= {u[0]}}, For[j = 1, j <= k, j = j + 1, l = {0}; 
For[i = 2, i <= n, i = i + 1, AppendTo[l, previous[dt, dx, {ans[[j]][[i - 1]], 
ans[[j]][[i]], ans[[j]][[i + 1]]}, alpha]]]; AppendTo[l, 0]; 
AppendTo[ans,l]]; ans[[k + 1]]]

However, it takes a lot of time to compute u[steps]=u[1000].And when want to plot the answer it gets stuck...

GraphicsGrid[Table[{Graphics[Table[{RGBColor[u[k][[i]], 0, 1 - u[k][[i]]], 
Rectangle[{points[[i]] - dx/2, 0}, {points[[i]] + dx/2, 1}]}, {i, 1, n}]]}, 
{k, 0, steps, IntegerPart[1/dt]}]]

What is the problem? Is there a better way to define u[k]?

Thank you for your help.

share|improve this question
Are you not using NDSolve on purpose ? –  b.gatessucks Mar 27 '14 at 14:54
Yes, my idea is to approximate the solution by numerical methods. –  Triana Mar 27 '14 at 14:58
But that's precisely what NDSolve does: approximate the solution by numerical methods. –  murray Mar 27 '14 at 15:05
I know, but I have to program the method, I can't use that function –  Triana Mar 27 '14 at 15:12
Could you edit this code can be cut and pasted into a Mathematica notebook and evaluated? Right now there are syntax errors and some definitions just aren't expressed as (InputForm) Mathematica. –  Pillsy Mar 27 '14 at 15:12

1 Answer 1

up vote 1 down vote accepted

OK, I think I can give you some tips about performance here. There are a couple things you do that really tend to slow you down, and which I would describe as Mathematica "anti-patterns". In particular, building arrays by repeatedly calling AppendTo is likely to be really slow (the time taken will grow quadratically in the length of the list), and accessing arrays in a loop by indexing into them also tends to be slow. You want to construct arrays in one go. For instance, though it's not a performance bottleneck, you can create your initial condition like so:

u0 = ReplacePart[1/(1 + points^2)^0.25, {1 -> 0., -1 -> 0.}];

Most arithmetic operations will work automatically with lists, element by element, so

1 + {1, 2, 3}
{2, 3, 4}  

This is going to be much faster than using a loop to do things. ReplacePart is there to zero out the first and last elements to match the boundary conditions once the rest is done.

Similarly, we can use the built-in function ListCorrelate to replace your function previous like so (I called mine next because it appears that you're time-stepping forwards):

next[dt_, dx_, alpha_, u0_] :=
  ListCorrelate[{0., 1., 0.} + alpha*(dt/dx^2)*{1, -2, 1}, u0], 1, 0.]

No loops, AppendTo or element-by-element access.

Finally, I generate the solution by just repeatedly applying next to the initial conditions and keeping a list of results, using the Mathematica function NestList. While I was at it, I decided to keep the points for all intermediate time steps, which your function u didn't seem to do (was that intentional)?

solve[dx_, dt_, alpha_, init_, steps_] :=
  NestList[next[dt, dx, alpha, #] &, init, steps];

My computer takes about 2.5 seconds to evaluate u[1000] using your definition; it evaluates solve[dx, dt, alpha, u0, 1000] in about 0.015 seconds. By using Mathematica's built in functions for list construction and transformation I was able to speed things up by more than a factor of 100, and my answer (by taking the last row of the result) is very close to yours (though it differs in the last couple decimal places due to the vagaries of floating point arithmetic).

Please note that I didn't attempt to check whether you've implemented your finite difference scheme correctly. It looks right at a glance, but I really recommend that you check your results against a set of initial conditions for which an exact solution exists.

share|improve this answer

Your Answer


By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.