OBLICZENIA CAŁEK FRESNELA

 PRZEZ ROZKŁAD NA

SZEREG FUNKCJI BESSELA

RZĘDU POŁÓWKOWEGO

 

 

MGR INŻ. JERZY KOSEK

 

 

 

Praca napisana na

Podyplomowym Studium Informatyki

Wydziału Budowy Maszyn i Informatyki

Akademia Techniczno – Humanistyczna

w Bielsku - Białej

 

 

 

Bielsko - Biała, maj 2004 r.


1.  Definicja całek Fresnela.

Celem niniejszej pracy było utworzenie programu obliczającego całki Fresnela. Procedura licząca całki została wykorzystana przeze mnie w programie obliczającym rozkład natężenia światła za przesłoną. Fresnel Augustin Jean (1788 – 1827), francuski fizyk, odegrał główną rolę w ugruntowaniu poglądu na falową naturę światła przez stworzenie teorii dyfrakcji i interferencji fal świetlnych. Teoria ta pozwala obliczyć rozkład natężenia światła w dowolnej odległości za przesłoną dowolnego kształtu, gdy światło rozchodzi się z zadanego rodzaju źródła. Wyniki obliczeń teoretycznych znakomicie zgadzają się z pomiarami doświadczalnymi w rzeczywistych układach optycznych. W teorii pojawiają się całki, nazwane od ich twórcy całkami Fresnela:

 

C(x) =                                cosinus Fresnela

S(x)  =                               sinus Fresnela

 

2.      Własności całek Fresnela.

Utworzenie odpowiedniego algorytmu obliczającego całki wymagało znalezienia właściwych metod matematycznych. Znaną własnością całek jest to, że można je rozwinąć w szeregi nieskończone:

 

C(x) = J1/2(π/2*x2) + J5/2(π/2*x2) + J9/2(π/2*x2)  + J13/2(π/2*x2) + ...

S(x)  = J3/2(π/2*x2) + J7/2(π/2*x2) + J11/2(π/2*x2) + J15/2(π/2*x2) + ...

 

gdzie Jn/2(π/2*x2) są funkcjami Bessela rzędu połówkowego [1].

 

Oznaczmy y = π/2*x2

Pierwsze z funkcji Jn/2 wyrażają się przez funkcje elementarne następująco:

 

J1/2(y) =

J-1/2(y) =

J3/2(y) = - +

J-3/2(y) = - -

 

 

Wzór ogólny ma postać:

 

Jn/2(y) =         n = 1,2,...

 

Powyższy wzór nie nadaje się do obliczeń numerycznych. Trzeba skorzystać z ogólnych własności funkcji walcowych, do których należą funkcje Bessela.

 

3.    Własności funkcji walcowych.

Funkcje walcowe rzędu n – tego są rozwiązaniami równania różniczkowego postaci:

 

w’’  + 1/y*w + (1 – n2/y2)*w = 0

 

gdzie w ogólności: y – liczba zespolona,    n – liczba rzeczywista

 

Wśród rozwiązań tego równania znajdują się funkcje spełniające związki rekurencyjne [2]:

 

wn+1(y)=2n/y*wn(y)-wn-1(y)        

 

Powyższe równanie spełniają między innymi interesujące nas funkcje walcowe Bessela rzędu połówkowego Jn/2(y). Stad otrzymamy użyteczny dla obliczeń numerycznych związek:

 

Jn+1(y)=2n/y*Jn(y)-Jn-1(y)        , gdzie   n = . . .

 

Dla n=1/2 z wzoru tego mamy:

 

 J3/2(y)=1/y*J1/2(y)-J-1/2(y)  

 

Podobnie dla n=3/2,5/2, 7/2 .. mamy:

 

J5/2(y)=3/y*J3/2(y)-J1/2(y)     

J7/2(y)=5/y*J5/2(y)-J3/2(y)     

J9/2(y)=7/y*J7/2(y)-J5/2(y)      itd.

 

 

 

4.    Wyprowadzenie algorytmu numerycznego.

Oznaczenia:

y:= π/2*x2

a = J 1/2(y) =

b = J-1/2(y) =

c = J 3/2(y) = i/y*a – b , gdzie i  = 1

 

W pętli programowej kolejne wartości Jn+1(y) będą obliczane wg schematu:

 

c= i/y*a – b

i = i + 2

b:=a;

a:=c;

i:=i+2;

 

Aby obliczyć całkę sinus Fresnela, trzeba sumować wartości funkcji Bessela dla i = 1, 5, 9, ..., dlatego piszemy instrukcję:

 

If  i mod 4 = 1 then sinF:= sinF + c

 

Natomiast dla obliczenia całki cosinus Fresnela trzeba sumować wartości funkcji Bessela dla i = 3, 7, 11, ..., dlatego piszemy instrukcję:

 

If  i mod 4 = 3 then cosF:=cosF + c

 

Powyższe instrukcje należy umieścić w pętli liczącej kolejne wartości „c”.

 

5.    Funkcja kończąca rekurencję.

Aby znaleźć właściwy warunek kończący pętlę rekurencyjną trzeba zauważyć, że kolejne wyrazy szeregów C(x), S(x) nie maleją monotonicznie. Danie prostego warunku „stopu”  typu:

 

While c  >min  do { kolejne instrukcje pętli programowej}

jest niewłaściwe, bo dla niektórych wartości „x” pewne wyrazy Jn+1(y) leżeć mogą blisko miejsc zerowych danej funkcji Bessela – wtedy rekurencja skończyłaby się zbyt szybko.

Wykorzystałem więc tę własność funkcji Bessela, że y = 0 jest jedynym wspólnym miejscem zerowym[3]. Dlatego właściwy warunek kończenia rekurencji wymaga sprawdzenia, czy co najmniej dwa kolejno dodawane wyrazy są jeszcze większe od zadanej dokładności:

While (i/y*a   > min ) and ( b > min ) do { kolejne instrukcje pętli programowej}

 

6.    Testowanie programu liczącego całki Fresnela.

Dla osiągnięcia dużej dokładności obliczeń nadałem zmiennym a,b,c typ extended. Procedurę przetestowałem dla parametru x Î <-100,100>,

w siatce punktów zestawionej w tabeli:

 

X

0 - 1

1 - 5

5 - 10

10 - 25

25 - 30

30 - 100

Krok

0.01

0.001

0.0005

0.0002

0.0001

0.0001

 

gdzie Krok – krok próbkowania w danym zakresie

 

Dla wartości x Î <-0.1, 0.1>  zbieżność procedury numerycznej uzyskałem dla min = 0. 000 001 , a dla x <= -0.1 lub x >= 0.1 dla min = 0. 000 000 1

 

Dla wartości „x” całek leżących między węzłami siatki zastosowałem interpolację liniową.

 

1.    Instrukcja do programu.

 

1)    Wykresy funkcji całkowych Fresnela.

Program rysuje wykresy funkcji całkowych:

C(x) =         dla x Î <-100,100>                                             

S(x)  =         dla x Î <-100,100>        

korzystając z własności, że:        

C(-x) = - C(x)          S(-x) = - S(x)         

a)    Dla x®±¥       C(x) ®±0.5,     S(x) ®±0. Pięknie widać to w programie: gdy za pomocą klawiszy „Przejdź dalej” lub „Przejdź w tył” zwiększamy  lub zmniejszamy wartość „x”, to wykresy funkcji zbliżają się do asymptot y=±0.5.

b)    Przebieg wykresów funkcji dla większych wartości „x” zagęszcza się i maleje do ±0.5; pola „SkalaY” i „Przyrost U/ pixel” pozwalają zwiększyć skalę wykresu w pionie i zmniejszyć krok próbkowania funkcji, przez co wykres powiększa wymiary. Widać wyraźnie, że nawet dla dużych wartości „x” obie funkcje oscylują wokół asymptot.

Uwaga:

Pola edycyjne „SkalaY” i „Przyrost U/ pixel” są aktywne, gdy „x” > 0.

 

 

2)    Obliczanie wartości całek

Program umożliwia obliczenie wartości całek Fresnela: zadajemy wartość  w polu edycyjnym „Wartość U” i naciskamy przycisk „Licz całki”. W polu „Ilość iteracji” podawana jest ilość kroków iteracyjnych, wykonanych przy obliczaniu całek. Widać, że ilość iteracji rośnie, gdy „x” rośnie, np. dla x = 30 ilość iteracji wynosi  2957.

 

 

  


 

 

Rys. 1 Przykładowe wyglądy panelu programu liczącego całki Fresnela.

 



[1] G. A. Korn, T. M. Korn, Matematyka dla pracowników naukowych i inżynierów,  PWN Warszawa 1983, s.276

[2] Tamże, s. 314

[3] Tamże, s. 318.