0
$\begingroup$

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]]]}
$\endgroup$
4
$\begingroup$

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}*) 
$\endgroup$
  • $\begingroup$ Thanks a lot! This is indeed a better way of defining these functions. $\endgroup$ – Pengfei Dec 28 '18 at 22:45
  • $\begingroup$ You're welcome. I'm still wondering about the errors you got ... $\endgroup$ – Ulrich Neumann Dec 29 '18 at 7:47
  • $\begingroup$ 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. $\endgroup$ – Pengfei Dec 29 '18 at 15:57

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.