# Optimization for iterating and plotting a matrix product and its elements?

I have an iterative matrix product of this type

$$\begin{pmatrix}y_{n+1} \\ x_{n+1} \end{pmatrix} = \begin{pmatrix} A & B \\ C & D \end{pmatrix} \begin{pmatrix}y_{n} \\ x_{n} \end{pmatrix}$$

Since I didn't know how to do it, I explicitly wrote the product as {A*yn + B*xn, C*yn + D*xn} and evaluated it with the NestList construct

Manipulate[
ListPlot[Orbit[d, F, R, B, k, n, pt],
PlotStyle -> {Black, PointSize[0]}, Frame -> True,
FrameTicks -> True, Axes -> False, AspectRatio -> 1,
ImageSize -> {600, 600}, ImagePadding -> 25, PlotRange -> Full],
"map parameters",
{{z, -0.7, "x0"}, -3, 3, .1, Appearance -> "Open",
ImageSize -> Tiny, AnimationRate -> 1/100},
{{w, -2.2, "y0"}, -3, 3, .1, Appearance -> "Open",
ImageSize -> Tiny, AnimationRate -> 1/100},
{{d, .012, "d"}, 0, 1, .001, Appearance -> "Open",
ImageSize -> Tiny, AnimationRate -> 1/100},
{{F, 0.044, "E"}, 0, 1, .001, Appearance -> "Open",
ImageSize -> Tiny, AnimationRate -> 1/10},
{{R, 0.76, "R"}, 0, 1, .001, Appearance -> "Open",
ImageSize -> Tiny, AnimationRate -> 1/10},
{{B, 6.157, "B"}, 0, 10, .001, Appearance -> "Open",
ImageSize -> Tiny, AnimationRate -> 1/10},
{{k, 10.04, "C"}, 0, 50, .001, Appearance -> "Open",
ImageSize -> Tiny, AnimationRate -> 1/10},
Button[Style["randomize parameters", 9], {z, w, d, F, R, B, k} =
RandomReal /@ {{0, 20}, {0, 20}, {0, 1}, {0, 1}, {0, 1}, {0,
10}, {0, 50}}, ImageSize -> 125, Appearance -> "FramedPalette"],
"\n number of iterates n",
{{n, 100, ""}, {100, 200, 1000, 5000, 100000, 500000}},
{{pt, {z, w}, ".."}},
ControlPlacement -> Left, FrameMargins -> 0,
Initialization :> (Orbit =
Compile[{d, F, R, B, k, {n, _Integer}, {pt, _Real, 1}},
NestList[
Module[{t =
B - k/(1 + #1[[1]] #1[[1]] + #1[[2]] #1[[2]])}, {1/
4 (5 F +
2 #1[[1]] (-1 + 2 R Cos[t] +
9 d R Sin[t]) + \[Sqrt]((F - 2 #1[[1]])^2 +
8 R #1[[1]] Sin[
t] (9 d F -
6 d #1[[1]] + (1 + 27 d^2) R #1[[1]] Sin[t]))) +
1/(2 R #1[[1]]) Csc[
t] (-F + 2 #1[[1]] -
12 d R #1[[1]] Sin[
t] - \[Sqrt]((F - 2 #1[[1]])^2 +
8 R #1[[1]] Sin[
t] (9 d F -
6 d #1[[1]] + (1 + 27 d^2) R #1[[1]] Sin[t]))) #1[[
2]], -R #1[[1]] Sin[t] + (#1[[2]]/(4 #1[[1]])) (
F + 2 #1[[
1]] (-1 - 2 R Cos[t] +
9 d R Sin[t]) + \[Sqrt]((F - 2 #1[[1]])^2 +
8 R #1[[1]] Sin[
t] (9 d F -
6 d #1[[1]] + (1 + 27 d^2) R #1[[1]] Sin[t])))}] &,
pt, n]])]


Here is a sample image showing my code output

I would also like to store the value of each element at every step, since I also didn't know how to do it I made a different program that first makes all the iterations and then I evaluate each element.

    Orbit = NestList[
Module[{t = o[#1[[1]], #1[[2]]]}, {#1[[1]] A[#1[[1]],#1[[2]],t] + #1[[2]] B[#1[[1]],#1[[2]],t], #1[[1]] C[#1[[1]],#1[[2]],t] + #1[[2]] D[#1[[1]],#1[[2]],t]}] &,
pt, n]])]
A[#1[[1]],#1[[2]],t]/. Thread[{x, y} -> #] & /@ %


Since my function is composed by many trigonometric functions I would like to ensure that the computer is using as many digits as possible. Is this done by default by Mathematica or I must add something else?

I know that a more efficient method should exist, but I'm out of my depth here.

• Your posted code includes symbols o,A,C,D in the complied function, but they are not defined anywhere. You should avoid using C and D anyway as they are built-in functions. – Simon Woods Apr 6 '16 at 20:03
• With A,B,C,D im refering to the matrix elements shown in the first image, those are functions of the variables #[[1]],#[[2]], o. o is also another function that depends on those variables, and it appears as the argument of the trigonometric functions that compose each element [A,B,C,D]. Each matrix element is a very long function so I replaced them with those symbols to avoid text cluttering. – Buffalo Apr 6 '16 at 21:36
• Maybe you can use some simpler, shorter functions to make a working example. It's really hard to optimise code just by looking at it - people really need something they can paste into Mathematica and play with. BTW you can use Dot for matrix multiplication. – Simon Woods Apr 6 '16 at 21:50
• I edited my post to add the full code. – Buffalo Apr 6 '16 at 22:03
• The square rooted term in the function appears 3 times and is quite complicated. It will be faster if you compute it just once like you do with t. Other than that I can't see any obvious optimisations. – Simon Woods Apr 7 '16 at 20:35