3
$\begingroup$

I need to iterate several functions such as {3 x[[1]]^2 + x[[2]], x[[2]], x[[3]]}, {Exp[x[[4]]], x[[2]], x[[3]], x[[4]]},etc. The number of variables depends on the problem I need to solve.

Here is an example:

ClearAll[aux]; 
aux[f_, α_, n_] := 
  Module[{g},
    g = Quiet @ Table[If[k == 1, f - α Random[] x[[k]], x[[k]]], {k, 1, n}];
    g]

ClearAll[fun]; 
fun[x_] = aux[3 x[[2]], 2, 3]

ClearAll[simul]; 
simul[f_, length_: 1, ic_] := Module[{i}, NestList[f[#] &, ic, length]]

simul[fun, 2, RandomReal[1, 3]]

Despite the warning message

Part::partd: Part specification x[[2]] is longer than depth of object

which can be suppressed, the above sequence of functions works (if somebody has a better solution, please let me know).

Unfortunately when I try to put everything together under a function, it won't work. For instance

trajec[f_, α_, n_] := 
  Module[{fun, orb},
    fun[x_] = Quiet @ aux[f, α, n];
    orb = simul[fun, n, RandomReal[1, n]];
    Table[{i, orb[[i]]}, {i, 1, n}]]

trajec[3 x[[2]], 2, 3]

returns

{1, {0.561078, 0.450353, 0.243782}}, 
{2, {-0.198531 x[[1]] + 3 x[[2]], x[[2]], x[[3]]}}, 
{3, {-0.198531 x[[1]] + 3 x[[2]], x[[2]], x[[3]]}}}

and not actual values.

What am I missing?

$\endgroup$
1
  • $\begingroup$ You will only get numeric results when you define a List x that contains numeric values otherwise x will be treated as a symbol. $\endgroup$
    – Matariki
    May 2 '15 at 20:51
2
$\begingroup$

What you are missing, I believe, is sufficient experience of Mathematica's core language at the functional level that you experimenting with. I give you credit for making a good try at formulating your code in a functional way, but I'm afraid you gone somewhat wide of the mark.

I have put your code into a form that works and which I think preserves your intent. I am not at all sure I have got it right because I find your intent obscure, so I hope you tell me if I am any closer to the mark than you were.

aux[f_, α_, x_List] :=
  Table[If[k == 1, f[x] - α RandomReal[] x[[k]], x[[k]]], {k, 1, Length[x]}]
fun[x_List] := aux[3 #[[2]] &, 2, x]
simul[f_, ic_, length_: 1] := NestList[f, ic, length]

SeedRandom[42]; simul[fun, RandomReal[1, 3], 2]
{{0.425905, 0.391023, 0.347069}, 
 {0.786568, 0.391023, 0.347069}, 
 {0.298463, 0.391023, 0.347069}}

Update

I do not recommend putting everything under one function. Rather, I would generalize the restricted case shown above as follows.

aux[f_, α_, x_] := 
 Table[If[k == 1, f[x] - α RandomReal[] x[[k]], x[[k]]], {k, 1, Length[x]}]

fun[α_, x_] := aux[3 #[[2]] &, α, x]

I write traject so that the length of the random vectors, len is decoupled from the number on iterations, n.

traject[f : (_Function | _Symbol), α_?NumericQ, 
        len_Integer /; len > 0, n_Integer /; n > 0] :=
  Module[{vec = RandomReal[1, len], i = 0, next, result},
    result = {0, vec};
    Nest[(result = {result, {++i, next = f[α, #]}}; next) &, vec, n];
    Cases[result, {_Integer, _List}, ∞]]

Now results such as the following can be computed.

 SeedRandom[42]; traject[fun, 2, 3, 3]
{{0, {0.425905, 0.391023, 0.347069}}, 
 {1, {0.786568, 0.391023, 0.347069}}, 
 {2, {0.298463, 0.391023, 0.347069}}, 
 {3, {1.00046, 0.391023, 0.347069}}}

The above is what you asked for. But many variants are possible, including the obvious: varying the number of iterations, varying the length of the generated vectors, and varying the value of α.

 SeedRandom[42]; traject[fun, 2, 3, 1]
{{0, {0.425905, 0.391023, 0.347069}}, 
 {1, {0.786568, 0.391023, 0.347069}}}
 SeedRandom[42]; traject[fun, π, 4, 2]
{{0, {0.425905, 0.391023, 0.347069, 0.453741}}, 
 {1, {0.429179, 0.391023, 0.347069, 0.453741}}, 
 {2, {0.783181, 0.391023, 0.347069, 0.453741}}}

More interesting, perhaps, is applying a different function

fun2[α_, x_] := Join[{RandomReal[α] E^x[[4]]}, Rest[x]]
SeedRandom[42]; traject[fun2, π, 4, 2]
{{0, {0.425905, 0.391023, 0.347069, 0.453741}}, 
 {1, {2.7495, 0.391023, 0.347069, 0.453741}}, 
 {2, {1.43008, 0.391023, 0.347069, 0.453741}}}
$\endgroup$
6
  • $\begingroup$ @m_goldberger. Many thanks. Yes, you are right when you say that I don't have sufficient experience of Mathematica's core language. Please believe when I say that I am trying but it hasn't been easy. The code you sent works as needed. Many thanks. Finally, how to put everything under one function ? (Final part of my question). $\endgroup$
    – Ed Mendes
    May 2 '15 at 22:29
  • $\begingroup$ @EdMendes. I have made an update that I hope will satisfy the request you made in your comment. $\endgroup$
    – m_goldberg
    May 3 '15 at 0:46
  • $\begingroup$ @m_goldberger. Thank you ever so much. Please correct me if I am wrong. My intent of having everything under one function was to use it within Manipulate so that I could vary alpha and the function under fun. Do you think that can be done with the functions above? One last question. Could you so kind to explain (_Function | _Symbol), please? $\endgroup$
    – Ed Mendes
    May 3 '15 at 1:27
  • $\begingroup$ @EdMendes. The vertical bar (|) is Alternative. The pattern f : (_Function | _Symbol) restricts the argument f to have a head that is either Symbol or Function. Presumably, when it is a symbol, that symbol will name a defined function. The head Function identifies any pure function. Placing this restriction on the first argument isn't fool-proof, but does catch many common errors. $\endgroup$
    – m_goldberg
    May 3 '15 at 3:44
  • $\begingroup$ @EdMendes. I don't see any reason why traject can't be called in a Manipulate expression with α values supplied through a control such as an Animator, the default control type for Manipulate expressions that looks like a slider. Every function called in a Manipulate expression doesn't have to be define inside that expression; it only has to known to kernel at the time the Manipulate expression is evaluated. $\endgroup$
    – m_goldberg
    May 3 '15 at 3:59

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.