# Too long computation

I am trying to solve this iterative task:

Here the upper index refers to time, lower refers to coordinate. It is an approximation of this PDE problem:

So I have for example t=10, I fix c=dt/dx^2=1/4, and varying n - number of segments in x coordinate I have t*4*N^2 number of segments in time.

So in my code, I just go through all segments of x to get all segments of second layer of time.

n = 70;
dx = 1/n;
dt = 4/n^2;
f = Table [Sin[i dx]*Pi^2, {i, 0, n}];
u = Table[0 x + 0 y, {x, 1000}, {y, (n)^2}];

c = 1/4;
For[i = 1, i < 49, i++,

For[j = 2, j < n, j++,
Part[u, j, i + 1] =
Part[u, j, i] (1 - 2 c) +
c (Part[u, j + 1, i] + Part[u, j - 1, i]) + dt *f[[j]]
]]



The problem is that the code is too slow to compute. I have no idea why? As I see, this is the easiest code for this task, so maybe it will take time for computation, but not as much as this one. I'll be very happy to know, is there something wrong with this code?

• Use ConstantArray[0,{1000,n^2}] to speed up definition of u, and if you don't need actual exact results, use numeric, e.g., using f = N@Table[...] will cut time several orders of magnitude. Finally, this can probably be vectorized, eliminating the loops, but I don't have the time right now to evaluate that. – ciao May 13 at 17:20
• @ciao thank you for your answer, but the speed problem is not in defining u, and I do need results about what is my function at certain x at certain t – Yura Holubeu May 13 at 17:24
• Easy way to solve is use NDSolve – Mariusz Iwaniuk May 13 at 17:30
• @MariuszIwaniuk That is not what I need!!!!! :) – Yura Holubeu May 13 at 17:36
• Yura, the point I think ciao and Mariusz are trying to make is that modifying Parts of lists in a For loop is not "the Mathematica" way, and that is one reason why this is so slow. – Marius Ladegård Meyer May 13 at 18:10

EDIT:

Okay, I've had a little bit of time to look at the PDE as well as your iterative solution. I don't believe the code below is set up correctly as there are a few mistakes.

First, if you want $$c = 1/4$$, then it should be $$dt = 1/(4 n^2)$$. Second, there's a $$\pi$$ missing in the definition of f. Thirdly, the array should be of size $$1/dt \times 1/dx$$ or $$4n^2 \times n$$. I've chosen here to make time increase as we move down the rows. It makes the most sense to me, but if you want it to be different you could always reverse the indices.

n = 20;
dx = 1/n;
dt = 1/4 dx^2;
f = N@Sin[Subdivide[n - 1] Pi]*Pi^2;
u = ConstantArray[0., {1/dt, 1/dx}];
c = dt/dx^2;
For[i = 1, i < Length[u], i++,
u[[i + 1, 2 ;; -2]] =
u[[i, 2 ;; -2]] (1 - 2 c) + c (u[[i, 3 ;;]] + u[[i, ;; -3]]) +
dt f[[2 ;; -2]]
]
ListPlot3D[
u,
AxesLabel -> {"x", "t", "u(x,t)"},
PlotRange -> Full
]


Now if we want to know the values for a particular timestep (say 1000), we can just ask for u[[1000]]. Note that the x and t axes in the above graph are in "steps" rather than being from 0 to 1.

That being said, the PDE that you show is easily solved both numerically and analytically. So you could also solve it with NDSolve or DSolve.

sol = First@DSolve[{
D[u[x, t], t] == D[u[x, t], {x, 2}] + Pi^2 Sin[Pi x],
u[x, 0] == 0
},
u[x, t],
{x, t}
]
Plot3D[
u[x, t] /. sol,
{x, 0, 1},
{t, 0, 1},
AxesLabel -> {"x", "t", "u(x,t)"}
]


$$\left\{u(x,t)\to \left(1-e^{-\pi ^2 t}\right) \sin (\pi x)\right\}$$

Looks the same as the previous graph to me. This solution has the advantage of providing the analytical equation. We can gain some intuition just by looking at the equation. For example, I can see that as $$t \rightarrow \infty$$, we're left with just $$\sin(\pi x)$$. In fact, we can see that this happens relatively quickly.

If this were a more complicated function that did not have a closed form, I would still recommend using NDSolve over an iterative equation. NDSolve has many advanced methods built-in, I would be surprised if a simple iteration worked better. Even for modest $$n = 400$$, the array u quickly becomes unwieldy at $$640000 \times 400$$ elements. Plotting this would take forever. You would need to sample only parts of u like ListPlot3D[u[[;; ;; 100]]] just to even plot it in a reasonable time. So you end up calculating a lot of points that aren't very useful. NDSolve is also quite fast in this case, taking only 0.016 seconds. With the result from NDSolve or DSolve you can also immediately do other things like take the derivative w.r.t. x or t.

ORIGINAL:

Making only the first 2 changes suggested by @ciao (ConstantArray and N), the code takes about 0.025 seconds to run on my machine.

n = 70;
dx = 1/n;
dt = 4/n^2;
f = N@Table[Sin[i dx]*Pi^2, {i, 0, n}];
u = ConstantArray[0., {1000, n^2}];

c = 1/4;
For[i = 1, i < 49, i++,
For[j = 2, j < n, j++,
Part[u, j, i + 1] =
Part[u, j, i] (1 - 2 c) +
c (Part[u, j + 1, i] + Part[u, j - 1, i]) + dt*f[[j]]]]


If we vectorize the second For loop, I get another factor of 3 speedup. I would have to play around with it to figure out if it's possible to vectorize both For loops, but I don't have time at the moment.

n = 70;
dx = 1./n;
dt = 4./n^2;
f = Sin[Range[0, n] dx]*Pi^2;
u = ConstantArray[0., {1000, n^2}];
c = 1/4;
For[i = 1, i < 49, i++,
u[[2 ;; n - 1, i + 1]] = u[[2 ;; n - 1, i]] (1 - 2 c) +
c (u[[3 ;; n, i]] + u[[1 ;; n - 2, i]]) +
dt*f[[2 ;; n - 1]]
]

• I'd recommend changing f to f = Sin[Range[0, n] dx]*Pi^2; rather than the explicit 70. (+1) – MarcoB May 13 at 19:31
• @MarcoB Oops, good catch, thanks! – MassDefect May 13 at 19:39
• yes, thank you, it is so good that it is working! If it is not hard for you, can you advise a command by which I can take all values u[x] at certain t? I feel that this is a simple command with lists, but... I can't figure out what it is... With it, it will be possible to Plot the solution, because now I have no idea how a function looks like... – Yura Holubeu May 13 at 21:35
• @YuraHolubeu I don't believe my original iterative code was correct. I've added some better code, and also show how to use DSolve if you're at all interested. – MassDefect May 13 at 23:42
• @MassDefect your update is really wonderful! But can you please explain, what does f[[2 ;; -2]] mean? How to understand it? – Yura Holubeu May 14 at 0:47