Dopasować okrąg

POSTAWIENIE PROBLEMU

Wyznaczyć 'najlepiej dopasowany' okrąg w punkcie W = [x(t0), y(t0)] do krzywej zadanej parametrycznie: x = x(t), y = y(t).

Np. środkowy okrąg wygląda na 'najlepiej dopasowany' spośród trzech na rysunku dla paraboli w jej wierzchołku W

[Maple Plot]

> "Jak sprecyzować pojęcie 'najlepiej dopasowany okrąg'? Jak go znaleźć?":

Pomocnicze procedury

> restart; with(plots); with(student)

> okr := [r*cos(t)+S[1], r*sin(t)+S[2], t = 0 .. 2*Pi...
okr := [r*cos(t)+S[1], r*sin(t)+S[2], t = 0 .. 2*Pi...
okr := [r*cos(t)+S[1], r*sin(t)+S[2], t = 0 .. 2*Pi...
okr := [r*cos(t)+S[1], r*sin(t)+S[2], t = 0 .. 2*Pi...

> opcje := t0 = -d .. d, view = [-d .. d, -d .. d], f...
opcje := t0 = -d .. d, view = [-d .. d, -d .. d], f...
opcje := t0 = -d .. d, view = [-d .. d, -d .. d], f...
opcje := t0 = -d .. d, view = [-d .. d, -d .. d], f...

Uwaga: w MAPLE'u są takie funkcje jak: odległość, iloczyn skalarny; jednak poniżej definiujemy WŁASNE (dla wygody?).

> odl := proc (u, w) options operator, arrow; sqrt((u...

PRZYKŁAD

dla paraboli y=x^2 w wierzchołku W (dla młodszego licealisty)

> W := [0, 0]; A := [a, a^2]; B := [-a, (-a)^2] W "pobliżu" wierzch. W paraboli wybieramy symetrycznie punkty A,B

> S := [0, s] Dla trójkąta WAB środek S okręgu opisanego leży oczywiście na osi OY .

> s = sqrt(a^2+(s-a^2)^2) Warunek SW = SA (zapisany analitycznie) pozwoli wyznaczyć współrzędną s środka S.

> s := solve(s = sqrt(a^2+(s-a^2)^2),s) Rozwiązujemy więc to równanie (na papierze równie łatwo):

s := 1/2+1/2*a^2

Zatem widać, że gdy A,B są blisko W , czyli gdy a jest "malutkie", to wartość s jest prawie równa 1/2.
Stąd "najlepiej dopasowany" okrąg ma środek w
S = [0, 1/2 ] i promień 1/2 (uruchom poniższą animację).

> r:=s:display({plot(x^2,x=-1.3..1.3,color=blue),animate({okr,seq([t*a*j,t*a^2,t=0..1],j=[-1,1]),[a*(1-2*t),a^2,t=0..1]},a=-1..0,color=red)});

[Maple Plot]

>

POMYSŁ I..

Środki P =[ xp, yp ] okręgów opisanych na trójkątach ABC dla punktów A,B,C tej krzywej, bliskich punktowi W, "powinny" przybliżać środek S1 szukanego okręgu.

Uwaga: w MAPLE'u są takie funkcje jak: odległość, iloczyn skalarny; jednak poniżej definiujemy WŁASNE (dla wygody?).

> odl := proc (u, w) options operator, arrow; sqrt((u...

> x := 'x'; y := 'y'; ta := 'ta'; tb := 'tb'; tc := '... ("techniczne wyczyszczenie" zmiennych)

> W := [x(t0), y(t0)]; A := [x(ta), y(ta)]; B := [x(t...

W := [x(t0), y(t0)]

A := [x(ta), y(ta)]

B := [x(tb), y(tb)]

C := [x(tc), y(tc)]

P := [xp, yp]

> symAB := il((A+B)/2-P,B-A) = 0; symBC := il((B+C)/2... symetralne odcinków AB i BC

symAB := (-xp+1/2*x(tb)+1/2*x(ta))*(-x(ta)+x(tb))+(...

symBC := (-xp+1/2*x(tc)+1/2*x(tb))*(-x(tb)+x(tc))+(...

> roz := solve({symAB, symBC},{xp, yp}) przecięcie symetralnych

roz := {xp = 1/2*(y(ta)*y(tc)^2+y(tc)*x(tb)^2+y(ta)...
roz := {xp = 1/2*(y(ta)*y(tc)^2+y(tc)*x(tb)^2+y(ta)...
roz := {xp = 1/2*(y(ta)*y(tc)^2+y(tc)*x(tb)^2+y(ta)...
roz := {xp = 1/2*(y(ta)*y(tc)^2+y(tc)*x(tb)^2+y(ta)...
roz := {xp = 1/2*(y(ta)*y(tc)^2+y(tc)*x(tb)^2+y(ta)...
roz := {xp = 1/2*(y(ta)*y(tc)^2+y(tc)*x(tb)^2+y(ta)...

> assign(roz) przypisanie zmiennym xp, yp obliczanych wartości (już MAPLE 'tak ma')

> S1 := [limit(limit(limit(xp,ta = t0),tb = t0),tc = ...
szukane S1 jest granicą P=[xp,yp]

S1 := [(`@@`(D,2)(y)(t0)*x(t0)*D(x)(t0)-D(y)(t0)*x(...
S1 := [(`@@`(D,2)(y)(t0)*x(t0)*D(x)(t0)-D(y)(t0)*x(...

TO Trochę straszne; D to operator pochodnej -- częste w MAPLE'u .

> r1 :=simplify(odl(W,S1)); promień 'najlepiej dopasowanego okręgu' jest oczywiście równy odległości punktów W i S

r1 := sqrt((D(x)(t0)^2+D(y)(t0)^2)^3/(-`@@`(D,2)(y)...

> S := S1; r := r1

Zastosowanie.

Definiujemy krzywą (parametrycznie):

> x:=t->t;
y:=t->cos(t)-3/4;

x := proc (t) options operator, arrow; t end proc

y := proc (t) options operator, arrow; cos(t)-3/4 e...

> ('S') = simplify(S); ('r') = simplify(r) teraz środek S i promień r wyznaczonego okręgu zależą (tylko) od t0

S = [(cos(t0)*t0-2*sin(t0)+sin(t0)*cos(t0)^2)/cos(t...

r = sqrt(-(cos(t0)^2-2)^3/cos(t0)^2)

Np. dla t0 = pi mamy:

> ('S(pi)') = eval(S,t0 = Pi), ('r(pi)') = eval(r,t0 ...
('S(pi)') = eval(S,t0 = Pi), ('r(pi)') = eval(r,t0 ...
('S(pi)') = eval(S,t0 = Pi), ('r(pi)') = eval(r,t0 ...

S(pi) = [Pi, -3/4], r(pi) = 1

[Maple Plot]

Na koniec trochę obrazków (szczegóły MAPLE'a są trochę 'okrutne' -- może można je pominąć?):

> d := 4; display({plot([krzywa, srodki],color = [blu...
d := 4; display({plot([krzywa, srodki],color = [blu...

[Maple Plot]

Proszę oglądnąć w przypadku innych krzywych (wystarczy zmienić definicję x(t) i y(t) )

POMYSŁ II.

Punkty Q =[ xq, yq ] przecięcia normalnych (tj. prostopadłych do krzywej) wystawionych w punktach A,W tej krzywej, gdy A jest bliskich punktowi W , "powinny" przybliżać środek S2 szukanego okręgu.

> x := 'x'; y := 'y'; ta := 'ta'; t0 := 't0'; xq := '... (to "techniczne wyczyszczenie" zmiennych)

> W := [x(t0), y(t0)]; A := [x(ta), y(ta)]; Q := [xq,...

W := [x(t0), y(t0)]

A := [x(ta), y(ta)]

Q := [xq, yq]

> styA := [D(x)(ta), D(y)(ta)]; styW := [D(x)(t0), D(... wektory styczne do krzywej w A i W

> normA := il(styA,Q-A) = 0; normW := il(styW,Q-W) = ... normalne w A i B

normA := D(x)(ta)*(-x(ta)+xq)+D(y)(ta)*(-y(ta)+yq) ...

normW := D(x)(t0)*(-x(t0)+xq)+D(y)(t0)*(-y(t0)+yq) ...

> roz := solve({normA, normW},{xq, yq}) przecięcie normalnych

roz := {yq = -(-D(x)(t0)*x(t0)*D(x)(ta)+D(x)(t0)*D(...
roz := {yq = -(-D(x)(t0)*x(t0)*D(x)(ta)+D(x)(t0)*D(...
roz := {yq = -(-D(x)(t0)*x(t0)*D(x)(ta)+D(x)(t0)*D(...
roz := {yq = -(-D(x)(t0)*x(t0)*D(x)(ta)+D(x)(t0)*D(...

> assign(roz) przypisanie zmiennym xq, yq obliczanych wartości (już MAPLE 'tak ma')

> S2 := ['limit(xq,ta = t0)', 'limit(yq,ta = t0)'] szukane S2 jest granicą Q

S2 := [limit(xq,ta = t0), limit(yq,ta = t0)]

> r2:='odl(W,S2)'; promień 'najlepiej dopasowanego okręgu' jest równy odległości punktów W i S

r2 := odl(W,S2)

> S := S2; r := r2

Zastosowanie.

Definiujemy krzywą parametrycznie:

> x:=t->t;
y:=t->t^2/3-3;

x := proc (t) options operator, arrow; t end proc

y := proc (t) options operator, arrow; 1/3*t^2-3 en...

> S := eval(S2); r := eval(r2) teraz środek S i promień r wyznaczonego okręgu zależą (tylko) od t0

S := [-4/9*t0^3, -3/2+t0^2]

r := 1/18*sqrt(972*t0^2+432*t0^4+64*t0^6+729)

Np. dla t0 = 1 mamy:

> ('S(1)') = eval(S,t0 = 1), ('r(1)') = eval(r,t0 = 1...
('S(1)') = eval(S,t0 = 1), ('r(1)') = eval(r,t0 = 1...
('S(1)') = eval(S,t0 = 1), ('r(1)') = eval(r,t0 = 1...

S(1) = [-4/9, -1/2], r(1) = 1/18*sqrt(2197)

[Maple Plot]

Na koniec trochę obrazków (szczegóły MAPLE'a są trochę 'okrutne' -- może można je pominąć?):

> display({plot([krzywa, srodki],color = [blue, green...

[Maple Plot]

Proszę oglądnąć w przypadku innych krzywych (wystarczy zmienić definicję x(t) i y(t) )

POMYSŁ III.

Dla punktu M tej krzywej (bliskiegu punktowi W ), znajdujemy okrąg przechodzący przez W i M o środku N =[ xn, yn ] na normalnej (prostopadłej) do krzywej wystawionej w punkcie W. Punkt N "powinien" przybliżać środek S3 szukanego okręgu.

> x := 'x'; y := 'y'; xn := 'xn'; yn := 'yn'; tm := '... ("techniczne wyczyszczenie" używanych zmiennych)

> W := [x(t0), y(t0)]; M := [x(tm), y(tm)]; N := [xn,...

W := [x(t0), y(t0)]

M := [x(tm), y(tm)]

N := [xn, yn]

> styW := [D(x)(t0), D(y)(t0)] wektor styczny do krzywej w punkcie W

> nor := il(styW,W-N) = 0; sym := il((M+W)/2-N,M-W) =... normalna w W oraz symetralna odcinka MW

nor := D(x)(t0)*(-xn+x(t0))+D(y)(t0)*(-yn+y(t0)) = ...

sym := (-xn+1/2*x(t0)+1/2*x(tm))*(-x(t0)+x(tm))+(-y...

> roz := solve({nor, sym},{xn, yn}) ich przecięcie

roz := {xn = 1/2*(2*y(t0)*x(t0)*D(x)(t0)-2*D(x)(t0)...
roz := {xn = 1/2*(2*y(t0)*x(t0)*D(x)(t0)-2*D(x)(t0)...
roz := {xn = 1/2*(2*y(t0)*x(t0)*D(x)(t0)-2*D(x)(t0)...
roz := {xn = 1/2*(2*y(t0)*x(t0)*D(x)(t0)-2*D(x)(t0)...
roz := {xn = 1/2*(2*y(t0)*x(t0)*D(x)(t0)-2*D(x)(t0)...

> assign(roz) przypisanie zmiennym xn, yn obliczanych wartości (już MAPLE 'tak ma')

> S3 := [limit(xn,tm = t0), limit(yn,tm = t0)] szukane S3 jest granicą N

S3 := [(`@@`(D,2)(y)(t0)*x(t0)*D(x)(t0)-D(y)(t0)*x(...
S3 := [(`@@`(D,2)(y)(t0)*x(t0)*D(x)(t0)-D(y)(t0)*x(...

> r3 := simplify(odl(W,S3)) promień 'najlepiej dopasowanego okręgu' jest oczywiście równy odległości punktów W i S

r3 := sqrt((D(x)(t0)^2+D(y)(t0)^2)^3/(-`@@`(D,2)(y)...

> S := S3; r := r3

Zastosowanie.

Definiujemy krzywą parametrycznie:

> x:=t->t;
y:=t->t*(t-1)*(t+2)/3+1;

x := proc (t) options operator, arrow; t end proc

y := proc (t) options operator, arrow; 1/3*t*(t-1)*...

> ('S') = simplify(S); ('r') = simplify(r) teraz środek S i promień r okręgu zależą (tylko) od t0

S = [-1/18*(-15*t0^2+24*t0-26+27*t0^6+54*t0^5-18*t0...
S = [-1/18*(-15*t0^2+24*t0-26+27*t0^6+54*t0^5-18*t0...

r = 1/18*sqrt((13+9*t0^4+12*t0^3-8*t0^2-8*t0)^3/(3*...

Np. dla t0 = -1 mamy:

> ('S(-1)') = subs(t0 = -1,S), ('r(-1)') = subs(t0 = ...
('S(-1)') = subs(t0 = -1,S), ('r(-1)') = subs(t0 = ...
('S(-1)') = subs(t0 = -1,S), ('r(-1)') = subs(t0 = ...

S(-1) = [-23/18, 5/6], r(-1) = 1/162*sqrt(125)*sqrt...

[Maple Plot]

Na koniec trochę obrazków (szczegóły MAPLE'a są trochę 'okrutne')::

> display({plot([krzywa, srodki],color = [blue, green...

[Maple Plot]

Proszę oglądnąć w przypadku innych krzywych (wystarczy zmienić definicję x(t) i y(t) )

PODSUMOWANIE

Czy te pomysły dają TEN SAM efekt? Czy S1 = S2 = S3 ???

Dla konkretnej krzywej:

> x:=t->t;
y:=t->t^3-t^2+3/4;

x := proc (t) options operator, arrow; t end proc

y := proc (t) options operator, arrow; t^3-t^2+3/4 ...

> 'S1'=S1;'r1'=r1;
'S2'=S2;'r2'=r2;
'S3'=S3;'r3'=r3;

S1 = [((6*t0-2)*t0-3*t0^2+2*t0-(3*t0^2-2*t0)^3)/(6*...

r1 = sqrt((1+(3*t0^2-2*t0)^2)^3/(-6*t0+2)^2)

S2 = [-1/2*t0^2*(-3+36*t0^2-8*t0-54*t0^3+27*t0^4)/(...

r2 = sqrt((t0+1/2*t0^2*(-3+36*t0^2-8*t0-54*t0^3+27*...
r2 = sqrt((t0+1/2*t0^2*(-3+36*t0^2-8*t0-54*t0^3+27*...

S3 = [((6*t0-2)*t0-3*t0^2+2*t0-(3*t0^2-2*t0)^3)/(6*...

r3 = sqrt((1+(3*t0^2-2*t0)^2)^3/(-6*t0+2)^2)

Można 'zmusić' MAPLE'a by sprawdził za nas czy S1 = S2 = S3 oraz r1 = r2 = r3 :

> 'S1 - S2'= simplify( S1 - S2);
'S2 - S3'= evala( S2 - S3);
'r1 - r2'= simplify( r1 - r2);
'r2 - r3'= simplify( r2 - r3);

S1-S2 = [0, 0]

S2-S3 = [0, 0]

r1-r2 = 0

r2-r3 = 0

Dla innej krzywej mamy:

> x:=t->t;
y:=t->t^4-2;

x := proc (t) options operator, arrow; t end proc

y := proc (t) options operator, arrow; t^4-2 end pr...

> 'S1'=simplify(S1),'r1'=simplify(r1);
'S2'=simplify(S2),'r2'=simplify(r2);
'S3'=simplify(S3),'r3'=simplify(r3);

S1 = [-2/3*t0*(-1+8*t0^6), 1/12*(1+28*t0^6-24*t0^2)...

S2 = [-16/3*t0^7+2/3*t0, 1/12*(1+28*t0^6-24*t0^2)/t...

S3 = [-2/3*t0*(-1+8*t0^6), 1/12*(1+28*t0^6-24*t0^2)...

Można 'zmusić' MAPLE'a by sprawdził za nas czy S1 = S2 = S3 oraz r1 = r2 = r3 :

> 'S1 - S2'= simplify( S1 - S2);
'S2 - S3'= simplify( S2 - S3);
'r1 - r2'= simplify( r1 - r2);
'r2 - r3'= simplify( r2 - r3);

S1-S2 = [0, 0]

S2-S3 = [0, 0]

r1-r2 = 0

r2-r3 = 0

A ogólnie? (Najpierw 'wykasujemy' definicje krzywej)

> x:='x'; y:='y';

x := 'x'

y := 'y'

> 'S1'=simplify(S1),'r1'=simplify(r1);
'S3'=simplify(S3),'r3'=simplify(r3);

S1 = [(`@@`(D,2)(y)(t0)*x(t0)*D(x)(t0)-D(y)(t0)*x(t...
S1 = [(`@@`(D,2)(y)(t0)*x(t0)*D(x)(t0)-D(y)(t0)*x(t...
S1 = [(`@@`(D,2)(y)(t0)*x(t0)*D(x)(t0)-D(y)(t0)*x(t...

S3 = [(`@@`(D,2)(y)(t0)*x(t0)*D(x)(t0)-D(y)(t0)*x(t...
S3 = [(`@@`(D,2)(y)(t0)*x(t0)*D(x)(t0)-D(y)(t0)*x(t...
S3 = [(`@@`(D,2)(y)(t0)*x(t0)*D(x)(t0)-D(y)(t0)*x(t...

Słabo widać, więc "zmuśmy" MAPLE'a by uprościł różnice:

> 'r1 - r3' = simplify(r1-r3);
'S1 - S3' = simplify(S1-S3);

r1-r3 = 0

S1-S3 = [0, 0]

Zera są dowodem tożsamości: S1 = S3 i r1 = r3 !!! TO JEST DOWÓD !!!

(Pomijamy tu dowód S1 = S2 i r1 = r2 .) Oczywiście otrzymane wzory ogólne są poprawne przy pewnych założeniach, czego tu nie będziemy analizować.

Uwaga: w podręcznikach analizy matematycznej znaleziony tu "najlepiej dopasowany" okrąg nazywa się okręgiem krzywizny , jego środek -- środkiem krzywizny , promień -- promieniem krzywizny , a odwrotność promienia -- krzywizną krzywej w punkcie W.