Habs nur kurz mit Maple gemacht, schaut noch ziemlich hässlich aus... leider nicht mehr Zeit gehabt (hoffe Du blickst durch).
p[0]: Startpunkt
g[0]: Startrichtung
n: Rekursionen
Hab die Variablen jetzt in jedem Schritt neu deklariert, falls Du dir die Kurve rauslassen willst.
Im Prinzip spiegelt man in jedem Schritt den Richtungsvektor der Gerade entlang dem Normalenvektor der Ellipse und berechnet daraus die neue Gerade.
restart;
> with(plots):
> a:=1;b:=2;n:=20;
> p[0]:=[1/sqrt(2)*a,1/sqrt(2)*b];
> g[0]:=[1,2];
> P:={plot([a*cos(s),b*sin(s),s=0..2*Pi],scaling=constrained)}:
>
> for k from 1 to n do
#Schnittpunkt Gerade Ellipse bestimmen
> t[k-1]:=solve({(p[k-1][1]+t*g[k-1][1])^2/a^2+
> (p[k-1][2]+t*g[k-1][2])^2/b^2=1,t<>0},{t});
> t[k-1]:=subs(t[k-1],t);
> p[k]:=[p[k-1][1]+t[k-1]*g[k-1][1],p[k-1][2]+t[k-1]*g[k-1][2]];
#s[k-1] so bestimmen, dass die Parametrisierungskurve in s[k-1] durch
den Schnittpunkt läuft (braucht man für den Normalenvektor der
Ellipse)
> s[k-1]:=solve({a*cos(s)=p[k][1],b*sin(s)=p[k][2]},s);
> s[k-1]:=subs(s[k-1],s);
#Spiegeln der Gerade entlang dem Normalenvektor n der Ellipse im
Schnittpunkt (nach g[k] = g[k-1] - 2*<g[k-1],n>/<n,n>*n)
> sp:= 1/((b*cos(s[k-1]))^2+(a*sin(s[k-1]))^2)*
> (-b*cos(s[k-1])*g[k-1][1]-a*sin(s[k-1])*g[k-1][2]);
> g[k]:=[g[k-1][1]+2*sp*(b*cos(s[k-1])),g[k-1][2]+2*sp*(a*sin(s[k-1]))];
#Plot
> P:=op(P);
> Q:=plot([p[k-1][1]+s*g[k-1][1],p[k-1][2]+s*g[k-1][2],s=0..t[k-1]],scaling=constrained);
> P:={P,Q};
> end:
>
> display(P);
>
