0

I experienced some difficulty to have the program do what I wanted. I defined a conditional recursive algorithm to evaluate the next pair of coordinates (psi[n+1],phi[n+1]) from (psi[n], phi[n]):

$\displaystyle (\psi[n+1],\phi[n+1])=\begin{cases}(\psi[n]+2\phi[n], \phi[n]),& \text{ if } L<\psi[n]+2\phi[n]< 2\pi-L,\\ &\\ (\psi_2[\psi[n], \phi[n]],\phi_2[\psi[n], \phi[n]]),& \text{ elsewhere. }\end{cases}$

Here $\psi_2(u,v)$ and $\phi_2(u,v)$ are two functions. These functions work fine if I give the value directly. For example, if (psi[n], phi[n])=(u,v)=(5.14159,0.1) (this appears later), then the functions give (psi[n+1], phi[n+1])=(psi2[u,v],phi2[u,v])=(2.81768,1.95294).

However, the recursive program fails to execute at the first nontrivial condition (that is, when L < Mod[psi[n] + 2*phi[n], 2 Pi] < 2 Pi - L). Please let know if I can make the following description more clear. Thanks!

The functions:

R = 100;
b = 95;
r = 10;
T = ArcCos[(R^2 + b^2 - r^2)/(2*b*R)]; (*the angle Psi_A*)
L = ArcCos[(R^2 - r^2 - b^2)/(2*b*r)]; (*the angle psi_A*)
Phi[u_, v_] := ArcCos[(r*Cos[v] + b*Cos[u + v])/R]; (*pullback angle on Gamma_R*)
Psi[u_, v_] := u + v - Phi[u, v] - 2*Pi; (*pullback position*)
n[u_, v_] := Floor[(T - Psi[u, v])/(2*Phi[u, v])]; (*reflection numbers on Gamma_R*)
phi2[u_, v_] := ArcCos[(R*Cos[Phi[u, v]] - b*Cos[Psi[u, v] + (2*n[u, v] + 1)*Phi[u, v]])/r]; (*the returing angle on Gamma_r*)
psi2[u_, v_] := Psi[u, v] + (2*n[u, v] + 1)*Phi[u, v] + phi2[u, v];  (*the returning position*)

The conditional recursive part (essentially the structure I desribed in the title, just with two coordinates):

psi[0] = Pi;
phi[0] = 0.1;
Do[psi[n + 1] = 
   If[L < Mod[psi[n] + 2*phi[n], 2 Pi], 
    If[Mod[psi[n] + 2*phi[n], 2 Pi] < 2 Pi - L, psi[n] + 2*phi[n], 
     psi2[psi[n], phi[n]]], psi2[psi[n], phi[n]]];
  phi[n + 1] = 
   If[L < Mod[psi[n] + 2*phi[n], 2 Pi], 
    If[Mod[psi[n] + 2*phi[n], 2 Pi] < 2 Pi - L, phi[n], 
     phi2[psi[n], phi[n]]], phi2[psi[n], phi[n]]], {n, 0, 11}];
   Do[Print[{n, psi[n], phi[n]}], {n, 0, 12}]

The first nine or ten steps works fine. However, it fails to evaluate n=11 from n=10. The output ((I skipped the results here))

{0,\[Pi],0.1}
{1,3.34159,0.1}
...
{10,5.14159,0.1} 
{11,4.28843 +0.953167 (1+2 10[5.14159,0.1])+ArcCos[1/10 (57.9104 -95 Cos[4.28843 +0.953167 (1+2 10[5.14159,0.1])])],ArcCos[1/10 (57.9104 -95 Cos[4.28843 +0.953167 (1+2 10[5.14159,0.1])])]}
{12,If[ArcCos[35/76]<Mod[4.28843 +0.953167 (1+2 10[5.14159,0.1])+3 ArcCos[1/10 (57.9104 -95 Cos[4.28843 +0.953167 (1+2 10[5.14159,0.1])])],2 \[Pi]],If[Mod[psi[n]+2 phi[n],2 \[Pi]]<2 \[Pi]-L,psi[n]+2 phi[n],psi2[psi[n],phi[n]]],psi2[psi[n],phi[n]]],If[ArcCos[35/76]<Mod[4.28843 +0.953167 (1+2 10[5.14159,0.1])+3 ArcCos[1/10 (57.9104 -95 Cos[4.28843 +0.953167 (1+2 10[5.14159,0.1])])],2 \[Pi]],If[Mod[psi[n]+2 phi[n],2 \[Pi]]<2 \[Pi]-L,phi[n],phi2[psi[n],phi[n]]],phi2[psi[n],phi[n]]]}
Pengfei
  • 105
  • 4

1 Answers1

4

Your program might be clearer if you define the functions recursive without Do

Clear[psi, phi]
psi[0] = Pi;
phi[0] = 0.1;

psi[n_] :=psi[n] = If[L < Mod[psi[n - 1] + 2*phi[n - 1], 2 Pi] < 2 Pi - L,psi[n - 1] + 2*phi[n - 1],psi2[psi[n - 1], phi[n - 1]]]
phi[n_] :=phi[n] = If[L < Mod[psi[n - 1] + 2*phi[n -1], 2 Pi] < 2 Pi - L,phi[n - 1], phi2[psi[n - 1], phi[n- 1]]]

Table[psi[i], {i, 1, 12}]
(*{3.34159, 3.54159, 3.74159, 3.94159, 4.14159,4.34159, 4.54159, 4.74159, 4.94159, 5.14159,2.81768, 2.7619}*)
Table[phi[i], {i, 1, 12}]
(*{0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1,1.95294, 1.16887}*)
Ulrich Neumann
  • 53,729
  • 2
  • 23
  • 55
  • Thanks a lot! This is indeed a better way of defining these functions. – Pengfei Dec 28 '18 at 22:45
  • You're welcome. I'm still wondering about the errors you got ... – Ulrich Neumann Dec 29 '18 at 07:47
  • It is probably caused by the repeated uses of 'n': one for the n[u_,v_] function in the first part of the code, and one for the recursive index. You used 'i' for the recursive index and the code works fine. – Pengfei Dec 29 '18 at 15:57