# How could I rewrite my code in professional style?

I reduced my question to the lowest possible for having your ideas:

ix[0] = 1;

vay[0] = 1;

h = 0.01 (*step size*)

ix[k_] := vay[k - 1] + ix[k - 1];

vay[k_] := ix[k - 1] - vay[k - 1];

x = Sum[h^n*ix[n], {n, 0, 5}];
y = Sum[h^n*vay[n], {n, 0, 5} ];


I must put:

x = ix[0] ;
y = vay[0];


and repeat it several times and finally plot x, y by considering step h. by using For loop and AppendTo I have got result for a small system but for a huge codes it will take more than 2 hours. I would be honored to see your suggestions about my problem.

RSolve is used to solve recurrence erquations

Clear[ix, vay]

eqns = {
ix[k] == vay[k - 1] + ix[k - 1],
vay[k] == ix[k - 1] - vay[k - 1],
ix[0] == 1, vay[0] == 1};

soln = RSolve[eqns,
{ix[k], vay[k]}, k][[1]] //
FullSimplify;

ix[k_] = ix[k] /. soln

(*  2^((1/2)*(-3 + k))*(2 + Sqrt[2] +
(-1)^k*(-2 + Sqrt[2]))  *)

vay[k_] = vay[k] /. soln

(*  2^(-1 + k/2)*(1 + (-1)^k)  *)


Verifying that these functions satisfy the equations

eqns // Simplify

(*  {True, True, True, True}  *)

h = 1/100;

Clear[x, y]

x[nmax_] = Sum[h^n ix[n], {n, 0, nmax}] // FullSimplify

(*  (1/4999)*(5100 +
(-102 - 101*Sqrt[2] +
(-1)^(1 + nmax)*(-102 +
101*Sqrt[2]))/
(2^((3/2)*(1 + nmax))*25^nmax))  *)


The sequence converges rapidly to

Limit[x[n], n -> Infinity]

(*  5100/4999  *)

% // N

(*  1.0202  *)

y[nmax_] = Sum[h^n vay[n], {n, 0, nmax}] // FullSimplify

(*  (5000 + (-100 - Sqrt[2] +
(-1)^(1 + nmax)*(-100 +
Sqrt[2]))/
(2^((3/2)*(1 + nmax))*25^nmax))/
4999  *)


The sequence converges rapidly to

Limit[y[n], n -> Infinity]

(*  5000/4999  *)

% // N

(*  1.0002  *)


EDIT:

Apparently you intend for ix[0] and vay[0] to vary between iterations, Therefore set them to variables rather than the constants {1, 1}.

Clear[ix, vay]
eqns = {ix[k] == vay[k - 1] + ix[k - 1], vay[k] == ix[k - 1] - vay[k - 1],
ix[0] == ix0, vay[0] == vay0};

soln = RSolve[eqns, {ix[k], vay[k]}, k][[1]] // FullSimplify;

ix[k_] = ix[k] /. soln

(*  2^((1/2)*(-3 + k))*
(ix0 + Sqrt[2]*ix0 +
(-1)^k*((-1 + Sqrt[2])*ix0 -
vay0) + vay0)  *)

vay[k_] = vay[k] /. soln

(*  2^((1/2)*(-3 + k))*
((1 + (-1)^(1 + k))*ix0 +
(-1 + Sqrt[2] + (-1)^k*
(1 + Sqrt[2]))*vay0)  *)

eqns // Simplify

(*  {True, True, True, True}  *)

h = 1/100;
nmax = 5;
Clear[x, y]

x[ix0_, vay0_] = Sum[h^n ix[n], {n, 0, nmax}] // FullSimplify

(*  (25005001*(101*ix0 + vay0))/
2500000000  *)

y[ix0_, vay0_] = Sum[h^n vay[n], {n, 0, nmax}] // FullSimplify

(*  (25005001*(ix0 + 99*vay0))/
2500000000  *)


To use the {x, y} output pairs as the {ix0, vay0} input pairs for the next iteration use NestList

xyData = NestList[{x @@ #, y @@ #} &, {1, 1}, 20] // N

(*  {{1., 1.}, {1.0202, 1.0002}, {1.04062, 1.0006}, {1.06124, 1.0012}, {1.08208,
1.002}, {1.10314, 1.003}, {1.12443, 1.00421}, {1.14594, 1.00561}, {1.16769,
1.00721}, {1.18968, 1.00902}, {1.21191, 1.01103}, {1.23439,
1.01324}, {1.25711, 1.01565}, {1.2801, 1.01827}, {1.30334,
1.0211}, {1.32685, 1.02412}, {1.35063, 1.02736}, {1.37469,
1.03079}, {1.39902, 1.03444}, {1.42364, 1.03829}, {1.44855, 1.04236}} *)

ListPlot[Transpose[xyData],
PlotLegends -> {"x", "y"},
DataRange -> {0, Length[xyData] - 1}]


EDIT 2:

To find functions that align with the x and y sequences, use FindSequenceFunction

xs[n_] = xyData[[All, 1]] // Rest // FindSequenceFunction[#, n] & //
Simplify

(*  (1/(100 + Sqrt[2]))*(2^(-1 - 8*n)*
(2*7^(2 + n)*
((3572143*(100 - Sqrt[2]))/
9765625)^n - 99*Sqrt[2]*
((25005001*(100 - Sqrt[2]))/
9765625)^n +
((25005001*(100 + Sqrt[2]))/
9765625)^n*(102 +
101*Sqrt[2])))  *)

ys[n_] = xyData[[All, 2]] // Rest // FindSequenceFunction[#, n] & //
Simplify

(*  2^(-1 - 8*n)*(25005001/9765625)^n*
((100 - Sqrt[2])^n +
(100 + Sqrt[2])^n)  *)

Plot[{xs[n], ys[n]}, {n, 0, 20},
Epilog -> {
AbsolutePointSize[4],
Blue,
Point[Transpose[{Range[0, 20], xyData[[All, 1]]}]],
Red,
Point[Transpose[{Range[0, 20], xyData[[All, 2]]}]]},
PlotLegends -> {"x", "y"}]


• Thank you so much for your comprehensive answer. But how can I repeat this process for having the next x and y (for the next steps)? (In each new process I have to put small “x” calculated from previous step as an initial condition for new step, i.e. x = ix[0] and y = vay[0].
– Sima
Nov 24, 2016 at 17:33
• You defined x and y as a Sum. I have derived the closed form of these sums as a function of the maximum value of the iterator. x[0] is ix[0] and y[0] is vay[0]. Presumably, your steps just increment the arguments to x and y, To see each step use x /@ Range[0, 5] // N and y /@ Range[0, 5] // N Nov 24, 2016 at 17:54
• I have tried to guess how you want the iterations performed. See edit. Nov 24, 2016 at 22:44
• Fantastic, Thank you so much.
– Sima
Nov 25, 2016 at 7:13
• While Mathematica does not analytically solve the problem by using RSlove, may you please suggest any other method to solve the problem?any other command to achieve what I would like to?Thank you @BobHanlon
– Sima
Nov 28, 2016 at 9:42

I have two main suggestions, though I think I might not be quite clear on what you are trying to achieve. First, in your definitions of ix and vay, if you add an extra part to the right hand side where you use Set as well as SetDelayed you don't need to recalculate every previous step every time, like this:

ix[k_] := ix[k] = vay[k - 1] + ix[k - 1];

vay[k_] := vay[k] = ix[k - 1] - vay[k - 1];


Next, using functional expressions like Map, Apply and Nest can do calculations significantly faster than loops like For. For example, here is a simple demonstration of a For loop vs a Map:

function[input_] := Sin[input] + Cos[input];

AbsoluteTiming[
For[
list1 = {}; iterator1 = 1,
iterator1 <= 10000,
iterator1++,
AppendTo[list1, function[iterator1]]
]
]

(*{2.24134, Null}*)

AbsoluteTiming[
list2 = function /@ Range[10000];
]
(*{0.0176008, Null}*)


Note: /@ is Map, and AbsoluteTiming tells you how many seconds it takes to calculate the result.

You can see the map is 200x faster, as well as being easier to write!

• Thank you for your answer. About your first suggestion: I have to recalculate each step in every round, because as you can see, the initial values are changing in each iteration, depending on the value of x and y, and maybe this is the problem which makes the process take very long time
– Sima
Nov 24, 2016 at 14:25
• I don't see where you change ix[0] and vay[0] in the code you posted. It might be useful if you posted your for loop. Nov 24, 2016 at 15:11