Mathematica Double Sequence Plotting in 3D

Dear friends,

Before talking about my problem, I would like to say I can code some programs in Visusal Basic 6.0, Matlab and Maple.

I have coded a program in Maple, but I have just learned that Mathematica draws much more better in 3D when compared to Maple.

So , I need to code the same program in Mathematica.
May I please get some information about the following topics:

• How to write a 'for next' loop?
• How to load values in to a double sequence (or a matrix, i.e. x(m,n) for m=1,2,...,20 and n=1,2,...,10)?
• How to plot a double sequence on 3D (with the form [m,n,x(m,n)])?

Thanks.

Note. If it could help I can also write the Maple codes.

Hi. Don't use For-Next loops. Use Table constructs. I'll show you how to set up a For loop and then show how it's implemented via a Table construct. Same dif with a 2-D array. The code below calculates at equally spaced values and then generates a ListPlot3D of the values.

Code:

---

sum = 0;
(* set up For loop to sum first 10 digits *)
For[i = 1, i <= 10, i++, sum += i; ];
sum

(* do same with applying "Plus" operator to a Table construct *)
Plus @@ Table[i, {i, 1, 10}]

(* set up 2D array of values for Sin(xy) and use ListPlot3D to plot it *)
slist = Table[{x, y, Sin[x*y]},
    {x, -1, 1, 0.2}, {y, -1, 1, 0.2}];
ListPlot3D[Flatten[slist, 1]]

---

shawsend ty vm.

If possible I will request more help (I feel myself like a lazzy student). (Blush)
I dont think that I can do sth like the following with these codes. (Shake)

Code:

---

#Maple codes that iterate a partial difference equation up to 20×20
restart:
with(plots):
with(plottools):
#Loading initial values
for m from -1 to -1 do
 for n from 1 to 20 do
  x[m,n]:=1:
 end do:
end do:
#Loading initial values
for m from -1 to 21 do
 for n from -1 to 0 do
  x[m,n]:=1:
 end do:
end do:
#Loading values for coefficient
for m from 0 to 20 do
 for n from 0 to 20 do
  if m=n then
  p[m,n]:=1:
  else
  p[m,n]:=1/((m+1)*(n+1)):
  end if:
 end do:
end do:
#Calculating the double sequence for 20×20
for i from 1 to 210 do
 u:=ceil((sqrt(8*i+1)-1)/2):
 m:=u*(u+1)/2-i:
 n:=u-m:
 x[m,n]:=evalf(x[m,n-1]-x[m+1,n-1]-p[m,n-1]*x[m-1,n-1-1]):
end do:
#Putting the values in to point form in R³
P:=[seq(seq(`if`(is(x[m,n],'real'),point([m,n,x[m,n]],color=`if`(x[m,n]<0,green,blue)),NULL),m=-1..20),n=-1..20)]:
#Plotting the double sequence
display(P,scaling=unconstrained,orientation=[40,50],projection=1,axes=boxed,symbolsize=8,symbol=circle);

---

http://img516.imageshack.us/img516/2471/outputuj1.jpg
If third component is negative, then I plot it in green; otherwise, it is blue.

Hi barkpuz. Mathematica indexes arrays from 1 to n. I noticed you're using -1 and 0. You'd have to convert your code to use arrays beginning with index 1 if you want to use list (arrays) in mathematica. Pretty sure that's the case anyway.

Quote:

---

Originally Posted by shawsend

Hi barkpuz. Mathematica indexes arrays from 1 to n. I noticed you're using -1 and 0. You'd have to convert your code to use arrays beginning with index 1 if you want to use list (arrays) in mathematica. Pretty sure that's the case anyway.

---

Thanks shawsend, now I changed the indices.
The following code iterates the following partial delay difference equation

where with all initial values set to .
These types of equations appear in Oscillation and Stability theory.

Code:

---

restart:
with(plots):
with(plottools):
#Number of points to iterate
S:=20:
#delay of the first component
tau:=3:
#Delay of the second component
sigma:=2:
#Loading iinitial values
for m from 1 to tau do
 for n from sigma+1 to (S+sigma+1) do
  x[m,n]:=1:
 end do:
end do:
#Loading initial values
for m from 1 to (S+tau+2) do
 for n from 1 to (sigma+1) do
  x[m,n]:=1:
 end do:
end do:
#Iterating the double sequence
for i from 1 to (S*(S+1)/2) do
 #Estimating values for triangle path
 u:=ceil((sqrt(8*i+1)-1)/2):
 m:=u*(u+1)/2-i:
 n:=u-m:
 #Calculating the coefficient
 A:=`if`((m+tau+1)=(n+sigma),1,1/((m+tau+1)*(n+sigma))):
 x[(m+tau+1),(n+sigma+1)]:=evalf(x[(m+tau+1),(n+sigma)]-x[(m+tau+2),(n+sigma)]-A*x[(m+1),n]):
end do:
#importing the values from the sequence




to build the point structure
P:=[seq(seq(`if`((m+n)<=(S+tau+sigma+2),point([(m-tau-1),(n-sigma-1),x[m,n]],color=`if`(x[m,n]<0,green,blue)),NULL),m=1..(S+tau+1)),n=1..(S+sigma+1))]:
#Plotting the point table
display(P,scaling=unconstrained,orientation=[40,50],projection=1,axes=boxed,symbolsize=8,symbol=circle);

---

gives the following graphic
http://img46.imageshack.us/img46/5008/32083252cq3.jpg
I would like to explain the details of this code.
For triangle path please see the following table.
http://img329.imageshack.us/img329/3631/81789017aq4.jpg
i.e. in the first step we are at , in the second step we are at , in the third one , and so on...
As you see from the graphic, the table has no values in the lower triangle, therefore I need to use to check if we are still at the upper triangle of the table (this can also be checked in different ways, i.e. `if`(is(x[m,n],'real'), but I guess that the one I now give in the code is more suitable for Mathematica).
point([m-tau-1,n-sigma-1,x[m,n]],color=`if`(x[m,n]<0,green,blue)) is to build the point structure in with the specified coordinates and color with respect to its value.

Do you think it is easy to convert this code into Mathematica?
Ty.

Hi Barkpuz. I think your code can be implemented in Mathematica. I've coded much, much more complex numerical PDE problems. However, I don't have the time to study your code in detail. I've made some attempt to convert it. I get array over-runs though. I'm doing something wrong. Perhaps you can further improve it. I know, it looks messy. Try and work with it though as this will help you learn Mathematica. Here's some notes:

1. Use Ceiling[arg] for the ceiling function.

2. Arrays use double brackets as in x[[m,n]].

3. Points are constucted as Point[{x,y,z}]

4. To code a colored point list, need the set of points in french brackets with a color preceding the point structure as in:

plist={{Blue,Point[{x1,y1,z1}]},{Red,Point[{x,y,z}]},...}

Then can use Show[Graphics3D[plist]] to show them

5. I converted the code below to "raw input form" in order to cut and paste it. If you wish, you can cut and paste it into Mathematica, select the cell, then choose Cell/Convert To/ Standard Form to recover most of the math format.

6. In Mathematica, never use upper case letters for user variables; may conflict with standard definitions. So I changed A to a.

Code:

---

s = 20;
tau := 3
sigma := 2
x = Table[1, {m, 1, s + tau + 1}, {n, 1, s + sigma + 1}];
For[i = 1, i <= (1/2)*s*(s + 1), i++,
  u = Ceiling[(1/2)*(Sqrt[8*i + 1] - 1)];
  m = (1/2)*u*(u + 1) - i; n = u - m;
  a = If[m + tau + 1 == n + sigma, 1,
    1/((m + tau + 1)*(n + sigma)); ];
  x[[m + tau + 1,n + sigma + 1]] =
    x[[m + tau + 1,n + sigma]] - x[[m + tau + 2,n + sigma]] -
    a*x[[m + 1,n]]; ]
plist = Table[If[m + n <= s + tau + sigma + 2,
    If[x[[m,n]] <= 0, {Green, Point[{m - tau - 1,
        n - sigma - 1, x[[m,n]]}]},
    {Blue, Point[{m - tau - 1, n - sigma - 1, x[[m,n]]}]}]],
  {m, 1, s + tau + 1}, {n, 1, s + sigma + 1}]
Show[Graphics3D[plist, PlotRange -> {{-25, 25}, {-25, 25},
    {-5, 5}}]]

---

Quote:

---

Originally Posted by shawsend

Hi Barkpuz. I think your code can be implemented in Mathematica. I've coded much, much more complex numerical PDE problems. However, I don't have the time to study your code in detail. I've made some attempt to convert it. I get array over-runs though. I'm doing something wrong. Perhaps you can further improve it. I know, it looks messy. Try and work with it though as this will help you learn Mathematica. Here's some notes:

1. Use Ceiling[arg] for the ceiling function.

2. Arrays use double brackets as in x[[m,n]].

3. Points are constucted as Point[{x,y,z}]

4. To code a colored point list, need the set of points in french brackets with a color preceding the point structure as in:

plist={{Blue,Point[{x1,y1,z1}]},{Red,Point[{x,y,z}]},...}

Then can use Show[Graphics3D[plist]] to show them

5. I converted the code below to "raw input form" in order to cut and paste it. If you wish, you can cut and paste it into Mathematica, select the cell, then choose Cell/Convert To/ Standard Form to recover most of the math format.

6. In Mathematica, never use upper case letters for user variables; may conflict with standard definitions. So I changed A to a.

Code:

---

s = 20;
tau := 3
sigma := 2
x = Table[1, {m, 1, s + tau + 1}, {n, 1, s + sigma + 1}];
For[i = 1, i <= (1/2)*s*(s + 1), i++,
  u = Ceiling[(1/2)*(Sqrt[8*i + 1] - 1)];
  m = (1/2)*u*(u + 1) - i; n = u - m;
  a = If[m + tau + 1 == n + sigma, 1,
    1/((m + tau + 1)*(n + sigma)); ];
  x[[m + tau + 1,n + sigma + 1]] =
    x[[m + tau + 1,n + sigma]] - x[[m + tau + 2,n + sigma]] -
    a*x[[m + 1,n]]; ]
plist = Table[If[m + n <= s + tau + sigma + 2,
    If[x[[m,n]] <= 0, {Green, Point[{m - tau - 1,
        n - sigma - 1, x[[m,n]]}]},
    {Blue, Point[{m - tau - 1, n - sigma - 1, x[[m,n]]}]}]],
  {m, 1, s + tau + 1}, {n, 1, s + sigma + 1}]
Show[Graphics3D[plist, PlotRange -> {{-25, 25}, {-25, 25},
    {-5, 5}}]]

---

---

I guess this will help me to fix the problem.
I will work on this and understant the structure, I hope.
Thanks.

Thanks again, now I give the working one below.

Code:

---

s = 10;
tau := 3;
sigma := 2;
x = Table[1, {m, 1, s + tau + 2}, {n, 1, s + sigma + 1}];
For[i = 1, i <= (1/2)*s*(s + 1), i++,
 u = Ceiling[(1/2)*(Sqrt[8*i + 1] - 1)];
 m = (1/2)*u*(u + 1) - i;
 n = u - m;
 a = If[m + tau + 1 == n + sigma, 1, 1/((m + tau + 1)*(n + sigma))];
 x[[m + tau + 1, n + sigma + 1]] = x[[m + tau + 1, n + sigma]] - x[[m + tau + 2, n + sigma]] - a*x[[m + 1, n]];
 ]
plist = Table[ If[m + n <= s + tau + sigma + 2, If[x[[m, n]] <= 0, {Green, Point[{m - tau - 1, n - sigma - 1, x[[m, n]]}]}, {Blue, Point[{m - tau - 1, n - sigma - 1, x[[m, n]]}]}]], {m, 1, s + tau + 1}, {n, 1, s + sigma + 1}];
Show[Graphics3D[plist, PlotRange -> {{-tau, s}, {-sigma, s}, {-20, 20}}]]

---

The codes above gives the following graphic
http://img505.imageshack.us/img505/7420/mathetj9.jpg

That's close but whenever you get a red background like that, an error occurred. Click on the x at the upper right corner and I get "NULL is not a graphics 3D primitive". If you check the array for x, we're getting Nulls there. Somehow, the indexes into the array are not correct. Getting Nulls in some entries so cannot construct a correct graphics primitive for some entries, which in this case are Point constructs. I'll try and look at it a little more sometime today.

Edit: I checked the array via MatrixForm[x] and it looks like it has numbers in each entry although I'm not sure it is in the form you desire. The errors then must be from over-extending it's bounds by the calculations done on m and n.

Got an idea: Enter this code (after you create x) and use a full-screen for Mathematica:

MatrixForm[x]

That gives you the matrix. Check that to see if the matrix is what you expect.

Ok. I think we're making progress. The nulls arose from the incorrect usage of Table I used to construct the point list "plist". The If statement created an entry if but did not have a default if that condition was not met. I added a Point[{0,0,0}] for that case. Not sure that's what you want. Here's the updated code:

Code:

---

s = 10;
tau := 3;
sigma := 2;
x = Table[1, {m, 1, s + tau + 2}, {n, 1, s + sigma + 1}];
For[i = 1, i <= (1/2)*s*(s + 1), i++,
 u = Ceiling[(1/2)*(Sqrt[8*i + 1] - 1)];
 m = (1/2)*u*(u + 1) - i;
 n = u - m;
 a = If[m + tau + 1 == n + sigma, 1, 1/((m + tau + 1)*(n + sigma))];
 x[[m + tau + 1, n + sigma + 1]] =
  x[[m + tau + 1, n + sigma]] - x[[m + tau + 2, n + sigma]] -
  a*x[[m + 1, n]];]
plist = Table[
  If[m + n <= s + tau + sigma + 2,
    If[x[[m, n]] <= 0, {Green,
      Point[{m - tau - 1, n - sigma - 1, x[[m, n]]}]}, {Blue,
      Point[{m - tau - 1, n - sigma - 1, x[[m, n]]}]}],
    Point[{0, 0, 0}]], {m, 1, s + tau + 1}, {n, 1, s + sigma + 1}];
Show[Graphics3D[plist,
  PlotRange -> {{-tau, s}, {-sigma, s}, {-20, 20}}]]

---

You do know about the interactive 3D graphics capabilities in ver 6 right? just hold down the left mouse button over the plot to rotate it. No red background either.

With the following code, I have the following graphic, but it is getting smaller as x takes larges values.
Is it possible to scale the plot with respect to its third component?
And something more, is it possible to set the angle of view with the code (for instance in my mapple code it is [40,50])?

Code:

---

s = 12;
tau := 3;
sigma := 2;
x = Table[1, {m, 1, s + tau + 2}, {n, 1, s + sigma + 1}];
xmin = 1;
xmax = 1;
For[i = 1, i <= (1/2)*s*(s + 1), i++,
 u = Ceiling[(1/2)*(Sqrt[8*i + 1] - 1)];
 m = (1/2)*u*(u + 1) - i;
 n = u - m;
 a = If[m + tau + 1 == n + sigma, 1, 1/((m + tau + 1)*(n + sigma))];
 c = x[[m + tau + 1, n + sigma]] - x[[m + tau + 2, n + sigma]] -
  a*x[[m + 1, n]];
 x[[m + tau + 1, n + sigma + 1]] = c;
 xmin = If [c < xmin, c, xmin];
 xmax = If [c > xmax, c, xmax];
 ]
plist = Table[
  If[m + n <= s + tau + sigma + 2,
    If[x[[m, n]] <= 0, {Green,
      Point[{m - tau - 1, n - sigma - 1, x[[m, n]]}]}, {Blue,
      Point[{m - tau - 1, n - sigma - 1, x[[m, n]]}]}],
    Point[{0, 0, 0}]], {m, 1, s + tau + 1}, {n, 1, s + sigma + 1}];
Show[Graphics3D[plist,
  PlotRange -> {{-tau, s}, {-sigma, s}, {xmin, xmax}}]]

---

http://img521.imageshack.us/img521/527/56255514bh4.jpg

Use BoxRatios:

Code:

---

Show[Graphics3D[plist,
  PlotRange -> {{-tau, s}, {-sigma, s}, {xmin, xmax}}],
 BoxRatios -> {1, 1, 1}, Axes -> True]

---

Not sure if that's a ver 6 feature though. Also can use ViewPoint->{x,y,z} to get a view point if you don't have ver 6. Do a f1 search with the cursor over this command to get help (put cursor over the command and press f1 right?)

Also, here's a tip: If you're really serious about learning Mathematica, hook up with the gang at Drexel:

Math Forum Discussions - comp.soft-sys.math.mathematica

Search the archives for a few weeks, say hour a day. Study what everyone is doing, cut and paste their answers into your Mathematica and learn what the code is doing, start asking questions, in 6 weeks, you'll be tops.(Wink)

Also, if you do decide to join that forum, always remember to convert your code to raw input form before posting it to the forum (you know, select code, choose Cell/convert to/ raw input form).