ゼルニケ多項式

ゼルニケ多項式は光学分野でよく使われています。

波面収差をゼルニケ多項式分解すると、各ゼルニケ多項式は独立した波面の形に対応しそれぞれが古典的な収差にも対応しており、収差成分をしることができます。

ゼルニケ (F. Zernike) はノーベル物理学賞を受賞した人で光波の位相差など波面 (等位相面) による解析をはじめた人です。ここではゼルニケ多項式についてMapleを使ってご紹介します。

[Maple Plot] [Maple Plot]

はじめに Maple の準備をしておきます。線形代数パッケージとプロット・パッケージを準備します。

>    restart:
with(plots):
with(linalg):

Warning, the name changecoords has been redefined

Warning, the protected names norm and trace have been redefined and unprotected

ゼルニケ多項式

ゼルニケ多項式は次のような複素関数で n と m の二つの次数を持ちます。

ただし実用的には実数関数として多く使用されています。

また単位円 (半径が 1) 上の関数で極座標の引数 ( rho  , theta  )をもちます。

すなわち範囲   0 <= rho _ <= 1 0 <= _ theta <= 2*Pi  です。

>    Zernike:=Z_r(n,m,rho)*exp(I*m*theta);

Zernike := Z_r(n,m,rho)*exp(m*theta*I)

ここで Z_r は次の様な関数です。

>    Z_r:=(n,m,rho)->sum( (-1)^k*((n-k)!/(k!*((n+m)/2-k)!*((n-m)/2-k)!))*rho^(n-2*k),
         k=0..1/2*(n-m) );

Z_r := proc (n, m, rho) options operator, arrow; sum((-1)^k*(n-k)!/k!/(1/2*n+1/2*m-k)!/(1/2*n-1/2*m-k)!*rho^(n-2*k),k = 0 .. 1/2*n-1/2*m) end proc

次数 m は -n から n まで2づつのステップでのみ定義されます。

実用上は実数関数として cos と sin の二つの式として使われます。

m が正のとき cos、負の時は sin 形式が使用されます。

m > 0 のとき (cos)

>    Z_cos:=(n,m,rho,theta)->   Z_r( n, abs(m), rho) * cos(m*theta);

Z_cos := proc (n, m, rho, theta) options operator, arrow; Z_r(n,abs(m),rho)*cos(m*theta) end proc

m < 0 のとき (sin)

>    Z_sin:=(n,m,rho,theta)-> Z_r( n, abs(m), rho) * sin(m*theta);

Z_sin := proc (n, m, rho, theta) options operator, arrow; Z_r(n,abs(m),rho)*sin(m*theta) end proc

                               cos         |             sin

n=0                               m=0

n=1                         m=1   m=-1

n=2                   m=2   m=0   m=-2

n=3                    m=3       m=1       m=-1    m=-3

n=4                m=4    m=2       m=0       m=-2    m=-4

次数 ( n = 2, m = 2 ) のプロットをしてみましょう。

>    plot3d( [r*cos(phi),r*sin(phi),Z_cos(2,2,r,phi)], r=0..1, phi=0..2*Pi,
   orientation=[75,70], axes=boxed, shading=ZHUE,grid=[50,50],  style=patchnogrid );
plot3d( [r*cos(phi),r*sin(phi),Z_cos(2,2,r,phi)], r=0..1, phi=0..2*Pi,
   orientation=[0,0], axes=none, shading=ZHUE,grid=[50,50], scaling=constrained,            style=patchnogrid,title= `Z_cos(2,2)`);

[Maple Plot]

[Maple Plot]

このように鞍点を持つ曲面となります。

>    plot3d( [r*cos(phi),r*sin(phi),Z_sin(2,2,r,phi)], r=0..1, phi=0..2*Pi,
     orientation=[75,70], axes=boxed, shading=ZHUE,grid=[50,50],   style=patchnogrid );
plot3d( [r*cos(phi),r*sin(phi),Z_sin(2,2,r,phi)], r=0..1, phi=0..2*Pi,
     orientation=[0,0], axes=none, shading=ZHUE,grid=[50,50],   scaling=constrained,            style=patchnogrid,title= `Z_sin(2,2)`);

[Maple Plot]

[Maple Plot]

このように同じ次数の場合、関数の形は同じで回転方向が45度傾いています。

光学におけるゼルニケ多項式

光学ではゼルニケ多項式を順番にならべ番号で扱われますその順番でゼルニケ多項式を生成してみましょう。

>    kk:=1:
for n from 0 to 8 do
for m from n by -2 to -n do
if m>=0 then
Z||kk:=Z_cos(n,m,r,phi);
else
Z||kk:=Z_sin(n,-m,r,phi);
fi;
print(kk,`      `,Z||kk);Z||kk:=unapply(Z||kk,r,phi);
kk:=kk+1:
od:od:

1, `      `, 1

2, `      `, r*cos(phi)

3, `      `, r*sin(phi)

4, `      `, r^2*cos(2*phi)

5, `      `, -1+2*r^2

6, `      `, r^2*sin(2*phi)

7, `      `, r^3*cos(3*phi)

8, `      `, (3*r^3-2*r)*cos(phi)

9, `      `, (3*r^3-2*r)*sin(phi)

10, `      `, r^3*sin(3*phi)

11, `      `, r^4*cos(4*phi)

12, `      `, (4*r^4-3*r^2)*cos(2*phi)

13, `      `, 1+6*r^4-6*r^2

14, `      `, (4*r^4-3*r^2)*sin(2*phi)

15, `      `, r^4*sin(4*phi)

16, `      `, r^5*cos(5*phi)

17, `      `, (5*r^5-4*r^3)*cos(3*phi)

18, `      `, (10*r^5-12*r^3+3*r)*cos(phi)

19, `      `, (10*r^5-12*r^3+3*r)*sin(phi)

20, `      `, (5*r^5-4*r^3)*sin(3*phi)

21, `      `, r^5*sin(5*phi)

22, `      `, r^6*cos(6*phi)

23, `      `, (6*r^6-5*r^4)*cos(4*phi)

24, `      `, (15*r^6-20*r^4+6*r^2)*cos(2*phi)

25, `      `, -1+20*r^6-30*r^4+12*r^2

26, `      `, (15*r^6-20*r^4+6*r^2)*sin(2*phi)

27, `      `, (6*r^6-5*r^4)*sin(4*phi)

28, `      `, r^6*sin(6*phi)

29, `      `, r^7*cos(7*phi)

30, `      `, (7*r^7-6*r^5)*cos(5*phi)

31, `      `, (21*r^7-30*r^5+10*r^3)*cos(3*phi)

32, `      `, (35*r^7-60*r^5+30*r^3-4*r)*cos(phi)

33, `      `, (35*r^7-60*r^5+30*r^3-4*r)*sin(phi)

34, `      `, (21*r^7-30*r^5+10*r^3)*sin(3*phi)

35, `      `, (7*r^7-6*r^5)*sin(5*phi)

36, `      `, r^7*sin(7*phi)

37, `      `, r^8*cos(8*phi)

38, `      `, (8*r^8-7*r^6)*cos(6*phi)

39, `      `, (28*r^8-42*r^6+15*r^4)*cos(4*phi)

40, `      `, (56*r^8-105*r^6+60*r^4-10*r^2)*cos(2*phi)

41, `      `, 1+70*r^8-140*r^6+90*r^4-20*r^2

42, `      `, (56*r^8-105*r^6+60*r^4-10*r^2)*sin(2*phi)

43, `      `, (28*r^8-42*r^6+15*r^4)*sin(4*phi)

44, `      `, (8*r^8-7*r^6)*sin(6*phi)

45, `      `, r^8*sin(8*phi)

これが最初の45個のゼルニケ多項式です。

いくつかのゼルニケ多項式をプロットしてみましょう。

>    print(Z4);
plot3d( [r*cos(phi),r*sin(phi),Z4(r,phi)], r=0..1, phi=0..2*Pi, orientation=[75,70], axes=boxed, shading=ZHUE,grid=[50,50],    style=patchnogrid );
plot3d( [r*cos(phi),r*sin(phi),Z4(r,phi)], r=0..1, phi=0..2*Pi, orientation=[0,0], axes=none, shading=ZHUE,grid=[50,50],   scaling=constrained,  style=patchnogrid,title= `Z4`);

proc (r, phi) options operator, arrow; r^2*cos(2*phi) end proc

[Maple Plot]

[Maple Plot]

>    print(Z5);
plot3d( [r*cos(phi),r*sin(phi),Z5(r,phi)], r=0..1, phi=0..2*Pi, orientation=[75,70], axes=boxed, shading=ZHUE,grid=[50,50],    style=patchnogrid );
plot3d( [r*cos(phi),r*sin(phi),Z5(r,phi)], r=0..1, phi=0..2*Pi, orientation=[0,0], axes=none, shading=ZHUE,grid=[50,50],   scaling=constrained,  style=patchnogrid,title= `Z5`);

proc (r, phi) options operator, arrow; -1+2*r^2 end proc

[Maple Plot]

[Maple Plot]

>        print(Z6);
plot3d( [r*cos(phi),r*sin(phi),Z6(r,phi)], r=0..1, phi=0..2*Pi, orientation=[75,70], axes=boxed, shading=ZHUE,grid=[50,50],    style=patchnogrid );
plot3d( [r*cos(phi),r*sin(phi),Z6(r,phi)], r=0..1, phi=0..2*Pi, orientation=[0,0], axes=none, shading=ZHUE,grid=[50,50],   scaling=constrained,  style=patchnogrid,title= `Z6`);

proc (r, phi) options operator, arrow; r^2*sin(2*phi) end proc

[Maple Plot]

[Maple Plot]

次は Z7 です。

>    print(Z7);
plot3d( [r*cos(phi),r*sin(phi),Z7(r,phi)], r=0..1, phi=0..2*Pi, orientation=[75,70], axes=boxed, shading=ZHUE,grid=[50,50],    style=patchnogrid );
plot3d( [r*cos(phi),r*sin(phi),Z7(r,phi)], r=0..1, phi=0..2*Pi, orientation=[0,0], axes=none, shading=ZHUE,grid=[50,50],   scaling=constrained,  style=patchnogrid,title= `Z7`);

proc (r, phi) options operator, arrow; r^3*cos(3*phi) end proc

[Maple Plot]

[Maple Plot]

Z8

>    print(Z8);
plot3d( [r*cos(phi),r*sin(phi),Z8(r,phi)], r=0..1, phi=0..2*Pi, orientation=[75,70], axes=boxed, shading=ZHUE,grid=[50,50],    style=patchnogrid );
plot3d( [r*cos(phi),r*sin(phi),Z8(r,phi)], r=0..1, phi=0..2*Pi, orientation=[0,0], axes=none, shading=ZHUE,grid=[50,50],   scaling=constrained,  style=patchnogrid,title= `Z8`);

proc (r, phi) options operator, arrow; (3*r^3-2*r)*cos(phi) end proc

[Maple Plot]

[Maple Plot]

ゼルニケ多項式の合成

ゼルニケ多項式を重み付けして足しあわすことで(線形結合)あらゆる曲面を生成することができます。

>    c2:=1:c10:=1:

>    f:=(r,phi)->c2*Z2(r,phi)+c10*Z10(r,phi);

f := proc (r, phi) options operator, arrow; c2*Z2(r,phi)+c10*Z10(r,phi) end proc

>    plot3d( [r*cos(phi),r*sin(phi),f(r,phi)], r=0..1, phi=0..2*Pi,
            orientation=[75,70], axes=boxed, shading=ZHUE,grid=[50,50],            style=patchnogrid );
plot3d( [r*cos(phi),r*sin(phi),f(r,phi)], r=0..1, phi=0..2*Pi,
            orientation=[0,0], axes=none, shading=ZHUE,grid=[50,50],            style=patchnogrid,title= `c2*Z2+c10*Z10`);

[Maple Plot]

[Maple Plot]

>    c2:=1:c15:=1:

>    f:=(r,phi)->c2*Z2(r,phi)+c15*Z15(r,phi);

f := proc (r, phi) options operator, arrow; c2*Z2(r,phi)+c15*Z15(r,phi) end proc

>    plot3d( [r*cos(phi),r*sin(phi),f(r,phi)], r=0..1, phi=0..2*Pi,
            orientation=[75,70], axes=boxed, shading=ZHUE,grid=[50,50],            style=patchnogrid );
plot3d( [r*cos(phi),r*sin(phi),f(r,phi)], r=0..1, phi=0..2*Pi,
            orientation=[0,0], axes=none, shading=ZHUE,grid=[50,50],            style=patchnogrid,title= `c2*Z2+c15*Z15`);

[Maple Plot]

[Maple Plot]

ゼルニケ多項式の直行性

ゼルニケ多項式は直行多項式系をなしておりこれが重要な役割を果たします。

すなわち自分自身以外のゼルニケ多項式とは直行しているので内積(単位円上での積分)がゼロとなります。

>    int(int(Z_cos(3,3,r,phi)*Z_cos(2,2,r,phi),phi=0..2*Pi)*r,r=0..1);

0

>    f:=(r,phi)->2*Z_cos(2,2,r,phi)-3*Z_cos(3,3,r,phi)+4*Z_cos(4,4,r,phi)-5*Z_cos(5,5,r,phi)+6*Z_cos(6,6,r,phi)-7*Z_cos(7,7,r,phi);

f := proc (r, phi) options operator, arrow; 2*Z_cos(2,2,r,phi)-3*Z_cos(3,3,r,phi)+4*Z_cos(4,4,r,phi)-5*Z_cos(5,5,r,phi)+6*Z_cos(6,6,r,phi)-7*Z_cos(7,7,r,phi) end proc

>    plot3d( [r*cos(phi),r*sin(phi),f(r,phi)], r=0..1, phi=0..2*Pi,
     orientation=[75,70], axes=boxed, shading=ZHUE,grid=[50,50],   style=patchnogrid );
plot3d( [r*cos(phi),r*sin(phi),f(r,phi)], r=0..1, phi=0..2*Pi,
       orientation=[0,0], axes=none, shading=ZHUE,grid=[50,50],   style=patchnogrid,title= `Z_cos(3,3,r,phi)-Z_cos(10,3,r,phi)`);

[Maple Plot]

[Maple Plot]

この関数とそれぞれのゼルニケ多項式との積分をとることで他のゼルニケ多項式との積分はゼロになるので

そのゼルニケ多項式だけの係数を取り出すことができます。ただし正規化するために係数をかける必要があります。

>    2*(2+1)/Pi*int(int(f(r,phi)*Z_cos(2,2,r,phi),phi=0..2*Pi)*r,r=0..1);

2

>    2*(3+1)/Pi*int(int(f(r,phi)*Z_cos(3,3,r,phi),phi=0..2*Pi)*r,r=0..1);

-3

>    2*(4+1)/Pi*int(int(f(r,phi)*Z_cos(4,4,r,phi),phi=0..2*Pi)*r,r=0..1);

4

>    2*(5+1)/Pi*int(int(f(r,phi)*Z_cos(5,5,r,phi),phi=0..2*Pi)*r,r=0..1);

-5

>    2*(6+1)/Pi*int(int(f(r,phi)*Z_cos(6,6,r,phi),phi=0..2*Pi)*r,r=0..1);

6

>    2*(7+1)/Pi*int(int(f(r,phi)*Z_cos(7,7,r,phi),phi=0..2*Pi)*r,r=0..1);

-7

このようにみごとに係数をみつけています。

しかし積分を精度よく計算することは難しいので実用性は少し劣ります。

最小二乗近似

ゼルニケ多項式による近似は最小二乗法によるほうが実用的です。線形結合なので

線形の近似で十分です。

簡単な例を示します。次のような線形結合を考えます。

>    f:=(r,phi)->a*Z5(r,phi)+b*Z6(r,phi)+c*Z7(r,phi);

f := proc (r, phi) options operator, arrow; a*Z5(r,phi)+b*Z6(r,phi)+c*Z7(r,phi) end proc

実験データとして既知の係数を決めておきます。

>    a:=1:b:=2:c:=3:   exact:=[[a],[b],[c]]:

簡単な例としてこのような単純な係数とします。

測定点に対応する変数のデータ・ポイントを用意します。

>    m:=5:n:=20:Lx:=vector(m):Lz:=vector(n):Lz:=matrix(m,n):

対応する関数値を求めます。

>    for j to m do
rr:=evalf(j/m):Lx[j]:=rr:
for k to n do
ph:=evalf(k*2*Pi/n):Ly[k]:=ph:
Lz[j,k]:=f(rr,ph):
od:od:

>    with(plottools):

Warning, the assigned name arrow now has a global binding

>    p:=NULL:

>    d:=0.05:

>    for j to m do
for k to n do
rr:=Lx[j]:ph:=Ly[k]:x:=evalf(rr*cos(ph)):y:=evalf(rr*sin(ph)):z:=evalf(f(rr,ph)):
p:=p,cuboid([x,y,z],[x+d,y+d,z+d/5]):
od:od:

データポイントを表示してみましょう。

>    display(p);

[Maple Plot]

すなわちこれが測定データあるいは解析対象データです。

近似するパラメータを未知のものとするため初期化します。

>    a:='a':b:='b':c:='c':

>    E:={seq(seq(f(Lx[j],Ly[k])=Lz[j,k],k=1..n),j=1..m)}:

決定したいパラメータのリスト{集合}を用意します。

>    V:= {a,b,c}:

これで準備はできましたので最小二乗近似を行います。

LINALG (線形代数) パッケージの leastsqrs を使用します。

>    ans:=leastsqrs(E, V);

ans := {a = .9999999993, b = 1.999999999, c = 3.000000002}

みごとに係数を見つけています (10 桁の精度)。

解の値を代入します。

>    assign(ans);

もとの(正確な)値と並べてみましょう。

>    evalm([[a],[b],[c]]=exact);

matrix([[.9999999993], [1.999999999], [3.000000002]]) = matrix([[1], [2], [3]])

もうひとつの近似例

もうひとつの例で数値実験してみましょう。

>    ff:=(r,phi)->a*Z11(r,phi)+b*Z12(r,phi)+c*Z13(r,phi)+d*Z14(r,phi)+e*Z15(r,phi)+f*Z16(r,phi);

ff := proc (r, phi) options operator, arrow; a*Z11(r,phi)+b*Z12(r,phi)+c*Z13(r,phi)+d*Z14(r,phi)+e*Z15(r,phi)+f*Z16(r,phi) end proc

実験データとして既知の係数を決めておきます。

>    a:=1:b:=-2:c:=3:d:=-4:e:=5:f:=-6:   exact:=[[a],[b],[c],[d],[e],[f]]:

簡単な例としてこのような単純な係数とします。

測定点に対応する変数のデータ・ポイントを用意します。

>    m:=5:n:=20:

対応する関数値を求めます。

>    for j to m do
rr:=evalf(j/m):Lx[j]:=rr:
for k to n do
ph:=evalf(k*2*Pi/n):Ly[k]:=ph:

Lz[j,k]:=ff(rr,ph):
od:od:

>    exact;

[[1], [-2], [3], [-4], [5], [-6]]

>    p:=NULL:

>    dd:=0.05:

>    for j to m do
for k to n do
rr:=Lx[j]:ph:=Ly[k]:x:=evalf(rr*cos(ph)):y:=evalf(rr*sin(ph)):z:=evalf(ff(rr,ph)):
p:=p,cuboid([x,y,z],[x+dd,y+dd,z+10*dd]):
od:od:

データポイントを表示しましょう。

>    display(p);

[Maple Plot]

すなわちこれが測定データあるいは解析対象データです。

近似するパラメータを未知のものとするため初期化します。

>    a:='a':b:='b':c:='c':d:='d':e:='e':f:='f':

>    E:={seq(seq(ff(Lx[j],Ly[k])=Lz[j,k],k=1..n),j=1..m)}:

決定したいパラメータのリスト{集合}を用意します。

>    V:= {a,b,c,d,e,f}:

これで準備はできましたので最小自乗近似を行います。

>    ans:=leastsqrs(E, V);

ans := {e = 4.999999990, a = 1.000000001, d = -3.999999998, c = 2.999999989, f = -6.000000001, b = -2.000000003}

値を代入します。

>    assign(ans);

元の値と近似値を比べてみましょう。

>    evalm([[a],[b],[c],[d],[e],[f]]=exact);

matrix([[1.000000001], [-2.000000003], [2.999999989], [-3.999999998], [4.999999990], [-6.000000001]]) = matrix([[1], [-2], [3], [-4], [5], [-6]])

>   

みごとに近似されています。

このように Maple の数式処理と数値計算がうまく活用できることがお分かりいただけると思います。