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 am a graduate student using a past student's Mathematica program. Will someone please explain to me exactly what this piece of code is doing? Specifically the Do[c[j] = {}, {j, 0, k}], Append command, and the Coefficient command. Thanks in advance.

inttstar[1] = (1/8)*(tstar[2] - tstar[0]);     (* Eq 2.31 *)

inttstar[i_] := inttstar[i] = (tstar[i + 1]/(i + 1) - tstar[i - 1]/(i - 1))/
    4 - (-1)^i*tstar[0]/(2 (i^2 - 1));    (* Eq 2.30 *)
inttstar[k];

Do[c[j] = {}, {j, 0, k}]
Do[Do[c[j] = Append[c[j], Coefficient[inttstar[j], tstar[i], 1]], {i, 
   0, k}], {j, 0, k}]

r = {};

Do[r = Append[r, c[i]], {i, 0, 15}]

r
share|improve this question
3  
Please be aware that your code was done by someone not used to program in Mathematica. You may find very ugly constructs if the code gets more complex afterwards –  belisarius Jan 22 at 1:51
    
It is building a list (in c) by appending a list of coefficients of tstart generated by intstar. Run it (assuming an appropriate k is set), then in a session evaluate intstar[somenumber], observe the result. Then evaluate Coefficient[%,tstar[somenumber]] and you'll see the coefficient associated with that (if it exists) tstar term. As an aside, this is not a good example of MM code, it would be a good exercise to re-write it. –  rasher Jan 22 at 1:54

1 Answer 1

up vote 10 down vote accepted

First off, it's apparent that k needs to take some (positive integer) value, since it's the end-value of the iterator later on. So I add:

k = 15;

The next bit of the code sets up a recursive function where inttstar[i] depends on periods i-1 and i+1. (This looks a lot like some economic model to be solved.) Notice the inttstar[i_]:= inttstar[i] = (* etc *) which ensures that existing results are saved and aren't recalculated every time.

inttstar[1] = (1/8)*(tstar[2] - tstar[0]);(*Eq 2.31*)

inttstar[i_] := 
  inttstar[i] = (tstar[i + 1]/(i + 1) - tstar[i - 1]/(i - 1))/4 - (-1)^i*
     tstar[0]/(2 (i^2 - 1));(*Eq 2.30*)

So for k=15, you get:

inttstar[k]

tstar[0]/448 + 1/4 (-(tstar[14]/14) + tstar[16]/16)

This is the bit that marks the original author out as someone who knows little about Mathematica programming. It's setting up some empty arrays to be filled in later.

Do[c[j] = {}, {j, 0, k}] 

Then they fill in the empty arrays with coefficients from the above equations.

Do[Do[c[j] = Append[c[j], Coefficient[inttstar[j], tstar[i], 1]], {i, 0, 
   k}], {j, 0, k}] 

For example:

c[14]

(* {-(1/390), 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -(1/52), 0, 1/60} *0

And then (for goodness sake!) they do the exact same thing of setting up another empty vector to store all the c[index] results. Notice the hard-coded 15 there, versus k further up. Very sloppy.

r = {};
Do[r = Append[r, c[i]], {i, 0, 15}]

The end result is a coefficient matrix

 r

{{1/2, 1/4, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, {-(1/8), 0, 1/8, 0, 0, 
  0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, {-(1/6), -(1/4), 0, 1/12, 0, 0, 0, 0, 0, 
  0, 0, 0, 0, 0, 0, 0},  (* etc *) 

So how can we improve this? The main thing to do would be to get rid of the Matlab/Fortran style song and dance routine with the empty lists that get appended to in Do loops. Notice the order of the iterators in the Table command. You then wouldn't have the c[j] variables defined, but judging by the result, all that was wanted was a coefficient matrix to go into some solution technique anyway. Yes, this one line can replace everything starting from that first Do loop.

Clear[r];
r = Table[Coefficient[inttstar[j], tstar[i], 1], {j, 0, k}, {i, 0, k}]

Some links you might find useful:

Avoiding procedural loops (common pitfalls)
Alternatives to procedural loops

The point is that it's not just an issue of speed and efficiency of code or what is more "Mathematica-ish" in coding style. This is an issue of readability. When you load your code up with all those loops and intermediate temporary variables, you make it harder to follow what the code is doing - exactly why you have found it difficult to understand what this past student's code does. Mathematica is my first and only programming language, so I find it hard to understand why people put up with all that hard-to-read procedural code.

As an aside, doing all this in exact numbers rather than floating point numbers will make the solution slower than you might like. I'd suggest adding an N[] command around the definition of r, so it would be r= N[Table[....(* etc *). I'm willing to bet you'll need to solve this numerically anyway using NSolve. This looks like a rational-expectations model in economics, and most of those don't have closed-form solutions.

share|improve this answer
1  
What a nice, patient, detailed answer! +1 –  rasher Jan 22 at 3:34
    
That is EVERY question I had about the code and then some! Thank you! –  gKirkland Jan 23 at 16:18
    
That is EVERY question I had about the code and then some! Thank you! The university I attend uses MATLAB in undergraduate courses, which is the reason this code looks MATLAB-ish. There are several blank arrays generated and then filled through the program, I attempted to replace them with the Table[] function as you have but I must have used the wrong syntax. You are correct, it is a program to solve some periodic non-linear dynamics problems (numerically). Once again thank you! You've saved me hours of trial-and-error code rewriting as this is my first experience with Mathematica. –  gKirkland Jan 23 at 16:24
    
Correct me if I'm wrong, I can replace: ` Do[c[j] = {}, {j, 1, k + 1}] Do[Do[c[j] = Append[c[j], Coefficient[pmult[j], tstar[i], 1]], {i, 0, k}], {j, 1, k}] s = {}; Do[s = Append[s, c[i]], {i, 1, k}] s = Transpose[s] ` With: ' s = Table[Coefficient[pmult[j], tstar[i], 1], {i, 0, k}, {j, 1, k}] s = Transpose[s] ' –  gKirkland Jan 23 at 16:37
    
@gKirkland yes in fact you can get rid of the s = Transpose[s] by swapping the iterators in the previous statement. –  Verbeia Jan 23 at 19:39

Your Answer

 
discard

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.