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'm trying to implement an upwind scheme for solving transport equation (at first, I work with no right hand side and the simpliest case possible - velocity is constant, initial data is smooth, etc). I have $$\partial_tf+\frac 1\varepsilon\partial_x f=0,\quad t>0,\quad x\in (0,1)$$ with initial and boundary data $$f(0,x)=f_{in}(x), \quad f(t,0)=0.$$

The classic explicit upwind scheme writes $$f(n+1,0)=0.$$ $$f(n+1,j)=f(n,j)\left(1-\frac{\Delta t}{\varepsilon \Delta x}\right)f(n,j)+\frac{\Delta t}{\varepsilon \Delta x}f(n,j-1)$$ With CFL (stability) condition $0\le \frac{\Delta t}{\varepsilon \Delta x}<1 $.

So, I write

HatFunction[x_] := If[Abs[x] < 1, Exp[1/(x^2 - 1)], 0];
InitialFunction[x_] := HatFunction[(x - 1/2)*6];(*initial data*)

n = 100;(*epsilon = 1/n; dx = 1/n, dt = 1/2n^2 , we want to have 2n^2 \
iterations*)
dx = 1/n;
eps = 1/n;
dt = 1/2/n^2;

InitData = ConstantArray[0, n + 1];
For[i = 1, i <= n + 1, i++, 
  InitData[[i]] = InitialFunction[(i - 1)/n]];

ccenter = ConstantArray[1 - dt/dx/eps, n + 1];
coffset = ConstantArray[dt/eps/dx, n + 1]; (*coefficients in the numerical scheme. *)
(*Here they're constant, but that will change in more complex models*)

Points = InitData;
For[i = 1, i <= n + 1, i++, Points[[i]] = (i - 1)/n]; (*coordinates of points in our mesh*)

Fp = InitData; (*Fp is the numerical solution*)


Niter = 10; (*number of iterations*)
Do[
  Fp = Fp*ccenter+RotateRight[Fp]*coffset;Fp[[1]]=0,      
{Niter}];
ListPlot[{Transpose[{Points, Fp}]}, Joined -> True]

The problem arises when I want to increase the number of iterations (after all, in order to get an approximate solution at time $T=1$ we want to make $2n^2$ iterations). If I set Niter=10, then everything works fine and the final plot is built instantly (didn't time it). However, if I set Niter=20 (or something bigger), then the evaluation wouldn't finish (I waited 5 minutes then aborted it), while the common sense suggests that the execution time should be linear with respect to the number of iterations.

Is there a way to improve the performance of this code? To optimize the list manipulation?

I'll be glad to hear all suggestions.

share|improve this question
1  
I see you are using numbers without any "dot" (e.g. 1 instead of 1.). This may force MMA to make exact computations and explode in computational complexity. You might try to add N to your expressions or a dot to any numeric constant –  user8074 Mar 30 at 19:19
    
@user8074 thanks for your advice, I'll give it a try tomorrow. –  TZakrevskiy Mar 30 at 21:14
    
@user8074 Indeed, your idea worked brilliantly! Everything works fine now. –  TZakrevskiy Apr 1 at 20:51

Your Answer

 
discard

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

Browse other questions tagged or ask your own question.