Use vpasolve to solve iterations of the same function when the number of iterations is N

Hi,first time asking a question here, I did search through but couldn't find an answer to this posted already.

I'm trying to cut down this code so instead of running it for each iteration of F(z), i can input the amount of iterations I want and it will compute just that.

You can see in the code that sols3 or sols6 are calculating F^N(z) - z for N = 3 and 6 respectively. I want to just use N to input this and not write code for sols3/4/5/6 etc.

The problem is that when N = 1 I want to vpasolve(F(z) - z, z)

but for N = 2 I want vpasolve(F(F(z)) - z, z)

N = 3, vpasolve(F(F(F(z))) - z, z)

if that makes sense.

I tried to make g = F^N(z) and vpasolve(g(z) - z, z)

but obviously thats wrong.

I need to do the same for the derivatives in the later part of the code too but I think if it's explained for the first part I should be able to do that part trivially.

Thank you very much :)

clear variables

%   Define symbolic variables r and phi.
%   Also need to define a numerical value for complex constant c...
syms z; r = 0.9; Omega = (sqrt(5)+1)/2; phi = 2*pi*Omega;
F = @(z) z^2 + r*exp(1i*phi)*z;
dFdz = @(z) 2*z + r*exp(1i*phi);

% To find N-periodic points
N = 7;
%g = F^N(z);


%   Obtain numerical approximations to the roots using symbolic algebra...
%   PERIOD-N POINTS (i.e., the fixed points of the map F^(N))

% solsN = vpasolve (g(z)- z,z);
% disp(['PERIOD-N POINTS:']); pN_points = vpa(solsN)
% size(pN_points,1)



%   PERIOD-3 POINTS [i.e., the fixed points of the map F^(3)]
sols3 = vpasolve(F(F(F(z))) - z,z);
disp(['PERIOD-3 POINTS:']); p3_points = vpa(sols3)
size(p3_points,1)

%   PERIOD-6 POINTS [i.e., the fixed points of the map F^(6)]
sols6 = vpasolve(F(F(F(F(F(F(z)))))) - z,z);
disp(['PERIOD-6 POINTS:']); p6_points = vpa(sols6)
size(p6_points,1)


%   Check for instability in a period-3 orbit...
disp([' ']); disp(['  CHECKING PERIOD-3 ORBIT STABILITY.'])
a3 = p3_points(3); F3_prime = dFdz(F(F(a3)))*dFdz(F(a3))*dFdz(a3);
fprintf("    For selected point: z = %s \n",string(a3));
fprintf("        |dF^(3)/dz| = %s \n",string(abs(F3_prime)));
if abs(F3_prime) > 1
    disp(['           ==> selected period-3 point is UNSTABLE!'])
elseif abs(F3_prime) < 1
    disp(['           ==> period-3 point is STABLE!'])
else
    disp(['           ==> period-3 point is MARGINALLY STABLE!'])
end

%   Check for instability in a period-6 orbit...
disp([' ']); disp(['  CHECKING PERIOD-6 ORBIT STABILITY.'])
a6 = p6_points(6); F6_prime = dFdz(F(F(F(F(F(a6))))))*dFdz(F(F(F(F(a6)))))*dFdz(F(F(F(a6))))*dFdz(F(F(a6)))*dFdz(F(a6))*dFdz(a6);
fprintf("    For selected point: z = %s \n",string(a6));
fprintf("        |dF^(6)/dz| = %s \n",string(abs(F6_prime)));
if abs(F6_prime) > 1
    disp(['           ==> selected period-6 point is UNSTABLE!'])
elseif abs(F6_prime) < 1
    disp(['           ==> period-6 point is STABLE!'])
else
    disp(['           ==> period-6 point is MARGINALLY STABLE!'])
end

Answer 1

[Voss]

One way is to define a function that takes F, N, and z as arguments and iterates N times to construct F(F(...F(z)...))-z. Then call that function with the appropriate value of N.

% function definition
function out = FN(F,N,z)
in = z;
for ii = 1:N
    in = F(in);
end
out = in - z;
end


% testing set-up
syms z;
r = 0.9;
Omega = (sqrt(5)+1)/2;
phi = 2*pi*Omega;
F = @(z) z^2 + r*exp(1i*phi)*z;

% compare original approach vs new approach with N = 3
sols3_original = vpasolve(F(F(F(z))) - z,z);
sols3_new      = vpasolve(FN(F,3,z),z);
isequal(sols3_original,sols3_new)

ans = logical
   1

% compare original approach vs new approach with N = 6
sols6_original = vpasolve(F(F(F(F(F(F(z)))))) - z,z);
sols6_new      = vpasolve(FN(F,6,z),z);
isequal(sols6_original,sols6_new)

ans = logical
   1

[Walter Roberson]

I was thinking that this would be a good use for fold but unfortunately fold() does not accept a symbolic function as the first parameter.

[Voss]

Looks like it has to be a function that accepts (at least) two inputs, so modifying F accordingly seems to work:

% testing set-up
syms z;
r = 0.9;
Omega = (sqrt(5)+1)/2;
phi = 2*pi*Omega;
F = @(z,~) z^2 + r*exp(1i*phi)*z;


% compare original approach vs new approach with N = 3
N = 3;
sols3_fold     = vpasolve(fold(F,z*ones(1,N+1))-z,z);
sols3_original = vpasolve(F(F(F(z))) - z,z);
isequal(sols3_fold,sols3_original)

ans = logical
   1

% compare original approach vs new approach with N = 6
N = 6;
sols6_fold     = vpasolve(fold(F,z*ones(1,N+1))-z,z);
sols6_original = vpasolve(F(F(F(F(F(F(z)))))) - z,z);
isequal(sols6_fold,sols6_original)

ans = logical
   1

[Walter Roberson]

Note that you are passing in a function handle to fold not a symbolic function. But still fold() gets the job done when used as you indicate.

[Voss]

You're correct. I was thinking you meant a function handle that accepts symbolic variables, rather than a symbolic function. My mistake. I've edited the terminology I used in my comment.