中心が原点ではなく，しかも，傾きのある２次元楕円は，５つのパラメータをもつ曲線

　　ｆ（ｘ，ｙ）＝ａ1ｘ^2 ＋ａ2ｘｙ＋ａ3ｙ^2＋ａ4ｘ＋ａ5ｙ−１＝０

の形に書けます．

　

　一般に，平面上の任意の位置にある５点が決まれば，ａ1〜ａ5が定まりますから，唯一の円錐曲線が決定されます．「円錐曲線は任意の５点で一意に定まる」は，パスカルの定理「円錐曲線，すなわち楕円，双曲線，放物線に内接する任意の六角形の三組の対辺の交点は同一直線上にある．」の重要な系になっています．→［補］

　

　しかし，ａ1〜ａ5が定まったとしても

　　ｆ（ｘ，ｙ）＝ａ1ｘ^2 ＋ａ2ｘｙ＋ａ3ｙ^2＋ａ4ｘ＋ａ5ｙ−１＝０

を馬鹿正直に描こうとすると，描画時間が大層長くなってしまうという欠点があります．

　

　そこで，楕円の基本情報，すなわち，長半径・短半径をそれぞれａ，ｂ，長軸とｘ軸のなす角度をθ，楕円の中心を（ｘ0，ｙ0）として，楕円モデル式を

　　｛（ｘ−ｘ0）ｃｏｓθ＋（ｙ−ｙ0）ｓｉｎθ｝^2／ａ^2＋｛（ｘ−ｘ0）ｓｉｎθ−（ｙ−ｙ0）ｃｏｓθ｝^2／ｂ^2＝１

とすると，楕円を媒介変数表示の形に表せるので，描画が短時間で済むようになります．

　

　今回のコラムでは，モデル式を最初から媒介変数表示しやすい形

　　ｆ（ｘ，ｙ）＝｛（ｘ−ｘ0）ｃｏｓθ＋（ｙ−ｙ0）ｓｉｎθ｝^2／ａ^2＋｛（ｘ−ｘ0）ｓｉｎθ−（ｙ−ｙ0）ｃｏｓθ｝^2／ｂ^2−１

として，多変量の非線形最小２乗法によるあてはめプログラムを紹介することにします．

　

＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝

　

［１］あてはめ例

　

＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝

　

［２］プログラムとその解説

　

　1460行では上記のモデル式

　　ｆ（ｘ，ｙ）＝｛（ｘ−ｘ0）ｃｏｓθ＋（ｙ−ｙ0）ｓｉｎθ｝^2／ａ^2＋｛（ｘ−ｘ0）ｓｉｎθ−（ｙ−ｙ0）ｃｏｓθ｝^2／ｂ^2−１

を定義していますが，このモデル式では未定係数が三角関数の中に含まれているため，どのように変形してもパラメータに関して線形とはなりません．すなわち，本質的に非線形なので，線形最小２乗法は使えないのです．したがって，多変量の非線形最小２乗法を用いる方法を採ることになります．→［３］参照

　

　また，1530行〜1670行では，１階偏微分

　　∂ｆ／∂ａ，∂ｆ／∂ｂ，∂ｆ／∂θ，∂ｆ／∂ｘ0，∂ｆ／∂ｙ0　　

を定義しています．差分近似ではなく，実際に偏微分した式を用いているので，差分に比べ，計算精度・計算速度ともに向上させることができます．

　

　非線形最小２乗法の計算では，まず最初に５個のパラメータａ，ｂ，θ，ｘ0，ｙ0の初期値が必要になります．ここで用いるサンプルデータ（1690行〜1950行）のように楕円の全周にわたって２元データ（ｘ，ｙ）が得られている場合は，かなり正確な初期値をプログラム上で設定することができます（2440行〜2660行）．

　

　この非線形モデルのあてはめに対しては，ガウス・ニュートン法を用いました（2570行〜4350行）．ガウス・ニュートン法では，計算が収束するまで数秒を要しましたが，あてはまりの良さに関しては満足のゆくものでした．

　

1000 '

1010 ' **** elliptic approximation ****

1020 '　　　　　　　　　　　　　2002/08/23 (C) ｻﾄｳ ｲｸﾛｳ

1030 SCREEN 3,0:CONSOLE ,,0,1

1040 CLS 3:WIDTH 80,25:COLOR 6

1050 GOSUB *INITIALIZE

1060 GOSUB *READ.DATA

1070 GOSUB *INIT.VAL

1080 GOSUB *EXECUTION

1090 GOSUB *SCATTER.DIAGRAM

1100 GOSUB *REPORT

1110 END

1120 '

1130 ' *** 初期設定 ***

1140 '

1150 *INITIALIZE:

1160 M=5:' No. of parameters

1170 DIM P(M),W(M)

1180 DIM P0(M),P1(M),DP(M)

1190 DIM AZ(M,M),BZ(M,M)

1200 DIM HESSE(M,M)

1210 DIM SD(M),SE(M)

1220 DIM COV(M,M)

1230 '

1240 NV=2:' No. of independent variables

1250 DIM X(NV)

1260 DIM AVE(NV),DEV(NV),SEM(NV)

1270 DIM XMAX(NV),XMIN(NV)

1280 FOR J=0 TO NV

1290　 XMAX(J)=-1E+10:XMIN(J)=1E+10

1300 NEXT J

1310 '

1320 RX=.7:RY=.8:BOX=1:KIND=1:POWER=0

1330 IX=.7:IY=.8:MARK=1:CLR=2:DRAW.COLOR=4

1340 TXT.COLOR=6:GRP.COLOR=5:REF.COLOR=6

1350 SCALE.KIND=1:PLOT=1

1360 PI=3.14159

1370 JX=.709:JY=1

1380 RX=IX*JX:RY=IY*JY

1390 '

1400 FORM$="(((X(1)-P(4))*COS(P(3))+(X(2)-P(5))*SIN(P(3)))/P(1))^2+(((X(1)-P(4))*SIN(P(3))-(X(2)-P(5))*COS(P(3)))/P(2))^2-1"

1410 RETURN

1420 '

1430 ' *** DEFINED EQUATION ***

1440 '

1450 *DEFINE.FORMULA:

1460 DEF FNY(X)=(((X(1)-P(4))*COS(P(3))+(X(2)-P(5))*SIN(P(3)))/P(1))^2+(((X(1)-P(4))*SIN(P(3))-(X(2)-P(5))*COS(P(3)))/P(2))^2-1

1470 RETURN

1480 '

1490 *DEFINE.FORMULA2:

1500 ON J GOSUB *J1,*J2,*J3,*J4,*J5

1510 RETURN

1520 '

1530 *J1:

1540 DEF FNJ(X)=-2*((X(1)-P(4))*COS(P(3))+(X(2)-P(5))*SIN(P(3)))^2/P(1)^3

1550 RETURN

1560 *J2:

1570 DEF FNJ(X)=-2*((X(1)-P(4))*SIN(P(3))-(X(2)-P(5))*COS(P(3)))^2/P(2)^3

1580 RETURN

1590 *J3:

1600 DEF FNJ(X)=2*((X(1)-P(4))*COS(P(3))+(X(2)-P(5))*SIN(P(3)))*(-(X(1)-P(4))*SIN(P(3))+(X(2)-P(5))*COS(P(3)))/P(1)^2+2*((X(1)-P(4))*SIN(P(3))-(X(2)-P(5))*COS(P(3)))*((X(1)-P(4))*COS(P(3))+(X(2)-P(5))*SIN(P(3)))/P(2)^2

1610 RETURN

1620 *J4:

1630 DEF FNJ(X)=-2*((X(1)-P(4))*COS(P(3))+(X(2)-P(5))*SIN(P(3)))*COS(P(3))/P(1)^2-2*((X(1)-P(4))*SIN(P(3))-(X(2)-P(5))*COS(P(3)))*SIN(P(3))/P(2)^2

1640 RETURN

1650 *J5:

1660 DEF FNJ(X)=-2*((X(1)-P(4))*COS(P(3))+(X(2)-P(5))*SIN(P(3)))*SIN(P(3))/P(1)^2+2*((X(1)-P(4))*SIN(P(3))-(X(2)-P(5))*COS(P(3)))*COS(P(3))/P(2)^2

1670 RETURN

1680 '

1690 *D:

1700 DATA 25

1710 DATA 2,0

1720 DATA 7,-2

1730 DATA 14,-3

1740 DATA 18,-4

1750 DATA 23,-4

1760 DATA 31,-3

1770 DATA 38,-1

1780 DATA 41,4

1790 DATA 40,8

1800 DATA 38,12

1810 DATA 34,13

1820 DATA 31,16

1830 DATA 26,17

1840 DATA 19,19

1850 DATA 13,20

1860 DATA 7,20

1870 DATA 2,19

1880 DATA -3,18

1890 DATA -6,17

1900 DATA -9,14

1910 DATA -10,12

1920 DATA -9,8

1930 DATA -7,5

1940 DATA -3,3

1950 DATA 0,1

1960 '

1970 ' *** データ読み込み ***

1980 '

1990 *READ.DATA:

2000 CLS 3

2010 RESTORE *D

2020 READ N1

2030 '

2040 DIM XX1(NV,N1),YH(N1),WT(N1)

2050 DIM Y1(N1),E1(N1)

2060 DIM JACOBI(M,N1)

2070 FOR J=1 TO N1:WT(J)=1:NEXT J

2080 '

2090 FOR I=1 TO N1

2100　 FOR J=1 TO NV

2110　　 READ XX1(J,I)

2120　 NEXT J

2130 NEXT I

2140 '

2150 PRINT "ただいまデータの取り込み中"

2160 FOR I=1 TO N1

2170 PRINT "*";

2180 Y1(I)=0

2190 X(0)=Y1(I)

2200 FOR J=1 TO NV

2210　 X(J)=XX1(J,I)

2220 NEXT J

2230 GOSUB *CALCULATE

2240 NEXT I

2250 '

2260 FOR J=0 TO NV

2270　 SS=DEV(J)-AVE(J)*AVE(J)/N1

2280　 DEV(J)=SQR(SS/(N1-1))

2290　 SEM(J)=DEV(J)/SQR(N1)

2300　 AVE(J)=AVE(J)/N1

2310 NEXT J

2320 RETURN

2330 '

2340 ' *** CALCULATE ***

2350 '

2360 *CALCULATE:

2370 FOR J=0 TO NV

2380　 IF XMAX(J)<X(J) THEN XMAX(J)=X(J)

2390　 IF XMIN(J)>X(J) THEN XMIN(J)=X(J)

2400　 AVE(J)=AVE(J)+X(J)

2410　 DEV(J)=DEV(J)+X(J)*X(J)

2420 NEXT J

2430 RETURN

2440 '

2450 ' ** 初期値 **

2460 '

2470 *INIT.VAL:

2480 X0=0:Y0=0

2490 FOR I=1 TO N1

2500　 X=XX1(1,I):Y=XX1(2,I)

2510　 X0=X0+X

2520　 Y0=Y0+Y

2530 NEXT I

2540 X0=X0/N1:Y0=Y0/N1

2550 P(4)=X0:P(5)=Y0

2560 '

2570 MAX=-1E+10:MIN=1E+10

2580 FOR I=1 TO N1

2590　 X=XX1(1,I):Y=XX1(2,I)

2600　 D=SQR((X-X0)^2+(Y-Y0)^2)

2610　 IF MAX<D THEN MAX=D:XMAX=X-X0:YMAX=Y-Y0

2620　 IF MIN>D THEN MIN=D

2630 NEXT I

2640 P(1)=MAX:P(2)=MIN

2650 IF XMAX=0 THEN P(3)=PI/2 ELSE P(3)=ATN(YMAX/XMAX)

2660 RETURN

2670 '

2680 ' *** 非線形最小２乗法 ***

2690 '

2700 *EXECUTION:

2710 PRINT:PRINT "ただいま計算中"

2720 GOSUB *GAUSS.NEWTON

2730 GOSUB *FINAL

2740 RETURN

2750 '

2760 ' *** GAUSS-NEWTON METHOD ***

2770 '

2780 *GAUSS.NEWTON:

2790 LMAX=50

2800 LOOP=0

2810 GOSUB *RESIDUAL:SO=SS

2820 GOSUB *TEMPORARY

2830 '

2840 GOSUB *JACOBIAN

2850 GOSUB *HESSIAN

2860 FOR I=1 TO M:P1(I)=P(I):NEXT I

2870 GOSUB *MARQUARDT

2880 GOSUB *UPDATE

2890 GOSUB *TRANSFORM

2900 '

2910 IF SN<SO THEN 2930

2920 FOR I=1 TO M:DP(I)=DP(I)/2:NEXT I:GOSUB *TRANSFORM:' [DAMPING]

2930 GOSUB *CONVERGENCE:IF SW=0 THEN RETURN

2940 SO=SN:IF LOOP<LMAX THEN 2840

2950 RETURN

2960 '

2970 ' ** X **

2980 '

2990 *JACOBIAN:

3000 FOR I=1 TO N1

3010　 FOR V=1 TO NV:X(V)=XX1(V,I):NEXT V

3020　 FOR J=1 TO M

3030　　 GOSUB *DEFINE.FORMULA2

3040　　 JACOBI(J,I)=FNJ(X)

3050　 NEXT J

3060 NEXT I

3070 GOSUB *JACOBIAN.SUB

3080 RETURN

3090 '

3100 ' ** Y **

3110 '

3120 *JACOBIAN.SUB:

3130 FOR I=1 TO N1

3140　 FOR V=1 TO NV:X(V)=XX1(V,I):NEXT V

3150　 GOSUB *DEFINE.FORMULA:YH=FNY(X)

3160　 E1(I)=Y1(I)-YH

3170　 YH(I)=YH

3180 NEXT I

3190 '

3200 ' ** XT*Y **

3210 '

3220 FOR J=1 TO M

3230　 S=0

3240　 FOR I=1 TO N1

3250　　 S=S+JACOBI(J,I)*E1(I)*WT(I)

3260　 NEXT I

3270　 W(J)=S

3280 NEXT J

3290 RETURN

3300 '

3310 ' ** XT*X **

3320 '

3330 *HESSIAN:

3340 FOR M1=1 TO M

3350　 FOR M2=M1 TO M

3360　　 S=0

3370　　 FOR I=1 TO N1

3380　　　 S=S+JACOBI(M1,I)*JACOBI(M2,I)*WT(I)

3390　　 NEXT I

3400　　 HESSE(M1,M2)=S:HESSE(M2,M1)=S

3410　 NEXT M2

3420 NEXT M1

3430 NZ=M

3440 RETURN

3450 '

3460 ' ** MARQUARDT **

3470 '

3480 *MARQUARDT:

3490 FOR I=1 TO M

3500　 FOR J=1 TO M:AZ(I,J)=HESSE(I,J):NEXT J

3510 NEXT I

3520 RETURN

3530 '

3540 ' ** UPDATE **

3550 '

3560 *UPDATE:

3570 GOSUB *INVERSE:' (XT*X)^(-1)

3580 GOSUB *UPDATE.SUB

3590 RETURN

3600 '

3610 ' ** (XT*X)^(-1)*XT*Y **

3620 '

3630 *UPDATE.SUB:

3640 FOR J=1 TO M

3650　 S=0

3660　 FOR K=1 TO M

3670　　 S=S+BZ(J,K)*W(K)

3680　 NEXT K

3690　 DP(J)=S

3700 NEXT J

3710 RETURN

3720 '

3730 ' ** TRANSFORM of PARAMETERS **

3740 '

3750 *TRANSFORM:

3760 FOR I=1 TO M

3770　 P(I)=P1(I)+DP(I)

3780 NEXT I

3790 GOSUB *RESIDUAL:SN=SS

3800 RETURN

3810 '

3820 ' ** SUM of RESIDUALS **

3830 '

3840 *RESIDUAL:

3850 SS=0

3860 FOR K=1 TO N1

3870　 FOR V=1 TO NV:X(V)=XX1(V,K):NEXT V

3880　 GOSUB *DEFINE.FORMULA:YH=FNY(X)

3890　 SS=SS+WT(K)*(Y1(K)-YH)^2

3900 NEXT K

3910 RETURN

3920 '

3930 ' *** CHECK of CONVERGENCE ***

3940 '

3950 *CONVERGENCE:

3960 'EPS=.0001

3970 EPS=.001

3980 SW=0

3990 IF ABS((SN-SO)/SO)>EPS THEN SW=1

4000 GOSUB *TEMPORARY

4010 RETURN

4020 '

4030 ' *** TEMPORARY OUTPUT ***

4040 '

4050 *TEMPORARY:

4060 LOOP=LOOP+1

4070 CLS 3

4080 PRINT "LOOP=";LOOP

4090 PRINT

4100 PRINT "PARAMETERS: "

4110 FOR I=1 TO M

4120　 PRINT "P(";I;")=";P(I)

4130 NEXT I

4140 PRINT

4150 PRINT "RESIDUALS=";SS

4160 RETURN

4170 '

4180 ' *** FINAL OUTPUT ***

4190 '

4200 *FINAL:

4210 CLS 3

4220 PRINT "LOOP=";LOOP

4230 PRINT

4240 PRINT "PARAMETERS: "

4250 FOR I=1 TO M

4260　 PRINT "P(";I;")=";P(I)

4270 NEXT I

4280 PRINT

4290 PRINT "RESIDUALS=";SN

4300 '

4310 GOSUB *PROPAGATION

4320 PRINT

4330 PRINT "何かキーを押してください"

4340 WHILE INKEY$="":WEND

4350 RETURN

4360 '

4370 ' ** INVERSE MATRIX ( AZ==>BZ ) **

4380 '

4390 *INVERSE:' [GAUSS-JORDAN]

4400 FOR I=1 TO NZ

4410　 FOR J=1 TO NZ:BZ(I,J)=0:NEXT J

4420　 BZ(I,I)=1

4430 NEXT I

4440 FOR I=1 TO NZ

4450　 IZ=I

4460　 IF AZ(I,I)=0 THEN 4470 ELSE 4530

4470　 IZ=IZ+1:IF IZ>NZ THEN RZ=1:RETURN

4480　 IF AZ(IZ,I)=0 THEN 4470

4490　 FOR J=1 TO NZ

4500　　 SWAP AZ(I,J),AZ(IZ,J)

4510　　 SWAP BZ(I,J),BZ(IZ,J)

4520　 NEXT J

4530　 AII=AZ(I,I)

4540　 FOR J=1 TO NZ

4550　　 AZ(I,J)=AZ(I,J)/AII

4560　　 BZ(I,J)=BZ(I,J)/AII

4570　 NEXT J

4580　 FOR K=1 TO NZ

4590　　 AKI=AZ(K,I):IF K=I THEN 4640

4600　　 FOR J=1 TO NZ

4610　　　 AZ(K,J)=AZ(K,J)-AKI*AZ(I,J)

4620　　　 BZ(K,J)=BZ(K,J)-AKI*BZ(I,J)

4630　　 NEXT J

4640　 NEXT K

4650 NEXT I

4660 RZ=0

4670 RETURN

4680 '

4690 ' *** 係数の誤差変換 ***

4700 '

4710 *PROPAGATION:

4720 FOR IS=1 TO M

4730　 SD(IS)=SQR(ABS(BZ(IS,IS)*SN/(N1-M)))

4740　 SE(IS)=SD(IS)/SQR(N1)

4750　 '

4760　 FOR JS=IS TO M

4770　　 COV(IS,JS)=BZ(IS,JS)*SN/(N1-M):COV(JS,IS)=COV(IS,JS)

4780　 NEXT JS

4790 NEXT IS

4800 RETURN

　

＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝

　

　4860行〜4960行では楕円が作図エリアのほぼ中央に描かれるように，オートスケーリングしています．これによって，見やすい位置に見やすい目盛り間隔でグラフを描いてくれるはずです．

　

　あてはめ例をみると，スケールの範囲はデータの範囲を偏りなくカバーしていますが，楕円の広がり具合がうまく表示されているとはいえません．スケール利用率の最大化のためには，4870〜4880行，4930〜4940行の方が好ましいのですが，その設定で描かれたグラフは窮屈そうにみえ，余裕をもたせるために，4890行・4950行のように設定しています．オートスケーリングの詳細については拙著「実験ノートのグラフ化技法と最新解析法」山海堂をご参照願います．

　

（問）２次元楕円

　　｛（ｘ−ｘ0）ｃｏｓθ＋（ｙ−ｙ0）ｓｉｎθ｝^2／ａ^2＋｛（ｘ−ｘ0）ｓｉｎθ−（ｙ−ｙ0）ｃｏｓθ｝^2／ｂ^2＝１

の変数ｘとｙにおいて，その取り得る最大値・最小値を求めよ．

　

4810 '

4820 ' *** SCATTER DIAGRAM ***

4830 '

4840 *SCATTER.DIAGRAM:

4850 CLS 3:WIDTH 80,25:COLOR 6

4860 'DMIN=XMIN(1):DMAX=XMAX(1)

4870 'PX=SQR(P(1)*P(1)*COS(P(3))*COS(P(3))+P(2)*P(2)*SIN(P(3))*SIN(P(3)))

4880 'DMIN=P(4)-PX:DMAX=P(4)+PX

4890 DMIN=P(4)-P(1):DMAX=P(4)+P(1)

4900 GOSUB *LINER.SCALE:UO=SMIN:UL=SMAX:DU=DIF:SU=STP

4910 '

4920 'DMIN=XMIN(2):DMAX=XMAX(2)

4930 'PY=SQR(P(1)*P(1)*SIN(P(3))*SIN(P(3))+P(2)*P(2)*COS(P(3))*COS(P(3)))

4940 'DMIN=P(5)-PY:DMAX=P(5)+PY

4950 DMIN=P(5)-P(1):DMAX=P(5)+P(1)

4960 GOSUB *LINER.SCALE:VO=SMIN:VL=SMAX:DV=DIF:SV=STP

4970 GOSUB *DRAW.AXIS

4980 GOSUB *DRAW.ELLIPSE

4990 GOSUB *TEMP

5000 '

5010 CLS 3

5020 RETURN

5030 '

5040 ' *** 楕円 ***

5050 '

5060 *DRAW.ELLIPSE:

5070 'DRAW.COLOR=4

5080 WUO=UO:WUL=UL

5090 WVO=VO:WVL=VL

5100 WINDOW(WUO,-WVL)-(WUL,-WVO)

5110 VIEW(SX1,SY1)-(SX2,SY2)

5120 '

5130 GOSUB *CONTOUR

5140 '

5150 WINDOW(0,0)-(639,399)

5160 VIEW(0,0)-(639,399)

5170 RETURN

5180 '

5190 ' *** 輪郭 ***

5200 '

5210 *CONTOUR:

5220 C=COS(P(3)):S=SIN(P(3))

5230 AAA=P(1)

5240 BBB=P(2)

5250 X0=P(4):Y0=P(5)

5260 '

5270 G=0

5280 FOR ANGLE=0 TO 2*PI+.01 STEP PI/180

5290　 GX=AAA*COS(ANGLE)*C-BBB*SIN(ANGLE)*S+X0

5300　 GY=AAA*COS(ANGLE)*S+BBB*SIN(ANGLE)*C+Y0

5310　 IF G=0 THEN PSET(GX,-GY),DRAW.COLOR:G=1

5320　 LINE-(GX,-GY),DRAW.COLOR

5330 NEXT ANGLE

5340 RETURN

5350 '

5360 ' *** 座標軸 ***

5370 '

5380 *DRAW.AXIS:

5390 CLS 3

5400 SX1=59:SY1=360*(1-RY):SX2=580*RX+59:SY2=360

5410 LINE(SX1,SY1)-(SX2,SY2),GRP.COLOR,B

5420 RR=3:GOSUB *X.SCALE2:GOSUB *X.REF2

5430 RR=3:GOSUB *Y.SCALE2:GOSUB *Y.REF2

5440 RR=3:GOSUB *DOT.PLOT

5450 GOSUB *TITLE.BACK2

5460 RETURN

5470 '

5480 ' *** 等分目盛り (X) ***

5490 '

5500 *X.SCALE2:

5510 CU=(SX2-SX1)/DU

5520 IF KIND=0 THEN KARA=-RR:MADE=0

5530 IF KIND=1 THEN KARA=0　:MADE=RR

5540 IF KIND=2 THEN KARA=-RR:MADE=RR

5550 FOR K=0 TO DU

5560　 LINE(CU*K+SX1,SY2+KARA)-(CU*K+SX1,SY2+MADE),GRP.COLOR

5570　 LINE(CU*K+SX1,SY1-KARA)-(CU*K+SX1,SY1-MADE),GRP.COLOR

5580 NEXT K

5590 RETURN

5600 '

5610 ' *** 参照値の記入 (X) ***

5620 '

5630 *X.REF2:

5640 REF.COLOR=6

5650 UW=UL-UO

5660 FOR K=0 TO DU STEP SU

5670　 F=UW*K/DU+UO:F$=STR$(F)

5680　 IF F>=0 THEN F$=MID$(F$,2)

5690　 FL=LEN(F$)

5700　 FOR L=1 TO FL

5710　　 GX=CU*K+59-FL*8/2+(L-1)*8:GY=23*16

5720　　 PUT(GX,GY),KANJI(ASC(MID$(F$,L,1))),PSET,REF.COLOR,0

5730　 NEXT L

5740 NEXT K

5750 RETURN

5760 '

5770 ' *** 等分目盛り (Y) ***

5780 '

5790 *Y.SCALE2:

5800 CV=(SY2-SY1)/DV

5810 IF KIND=0 THEN KARA=0　:MADE=RR

5820 IF KIND=1 THEN KARA=-RR:MADE=0

5830 IF KIND=2 THEN KARA=-RR:MADE=RR

5840 FOR K=0 TO DV

5850　 LINE(SX1+KARA,SY2-CV*K)-(SX1+MADE,SY2-CV*K),GRP.COLOR

5860　 LINE(SX2-KARA,SY2-CV*K)-(SX2-MADE,SY2-CV*K),GRP.COLOR

5870 NEXT K

5880 RETURN

5890 '

5900 ' *** 参照値の記入 (Y) ***

5910 '

5920 *Y.REF2:

5930 REF.COLOR=6

5940 VW=VL-VO

5950 FOR K=0 TO DV STEP SV

5960　 F=VW*K/DV+VO:F$=STR$(F)

5970　 IF F>=0 THEN F$=MID$(F$,2)

5980　 FL=LEN(F$)

5990　 FOR L=1 TO FL

6000　　 GX=48-FL*8+(L-1)*8:GY=360-CV*K-8

6010　　 PUT(GX,GY),KANJI(ASC(MID$(F$,L,1))),PSET,REF.COLOR,0

6020　 NEXT L

6030 NEXT K

6040 RETURN

6050 '

6060 ' *** プロット ***

6070 '

6080 *DOT.PLOT:

6090 WUO=UO:WUL=UL

6100 WVO=VO:WVL=VL

6110 WX=WUL-WUO

6120 BX=(WUL*SX1-WUO*SX2)/WX

6130 AX=(SX2-SX1)/WX

6140 WY=WVL-WVO

6150 BY=(WVL*SY2-WVO*SY1)/WY

6160 AY=(SY1-SY2)/WY

6170 '

6180 FOR I=1 TO N1

6190 X=XX1(1,I):Y=XX1(2,I)

6200 XX=AX*X+BX

6210 YY=AY*Y+BY

6220 CIRCLE(XX,YY),RR,1

6230 PAINT(XX,YY),CLR,1

6240 CIRCLE(XX,YY),RR,CLR

6250 NEXT I

6260 RETURN

6270 '

6280 ' *** タイトル ***

6290 '

6300 *TITLE.BACK2:

6310 LX=INT((580*RX+59)/8)

6320 LY=INT(360*(1-RY)/16)

6330 LZ=INT((290*RX+59)/8)

6340 L1$="elliptic approximation by "

6350 L2$="GAUSS-NEWTON METHOD"

6360 L$=L1$+L2$

6370 LOCATE LZ-LEN(L$)\2,LY-3:PRINT L$

6380 '

6390 LOCATE LX+3,LY:PRINT "PARAMETERS: "

6400 FOR K=1 TO M

6410　 V$="P"+MID$(STR$(K),2)+"="

6420　 LOCATE LX+3,LY+K+1:PRINT V$;:PRINT P(K)

6430 NEXT K

6440 RETURN

6450 '

6460 *TEMP:

6470 LOCATE 0,24:PRINT "何かキーを押してください";:WHILE INKEY$="":WEND

6480 RETURN

6490 '

6500 ' *** 平均律スケールの作成 ***

6510 '

6520 *LINER.SCALE:

6530 DIM INTERVAL(10)

6540 RESTORE *REF.INTERVAL

6550 FOR I=1 TO 10:READ INTERVAL(I):NEXT I

6560 MAX.INTERVAL=10

6570 '

6580 ARGUMENT=DMAX-DMIN

6590 GOSUB *SCALE.RTN

6600 INDEX=CMP

6610 INTERVAL=INTERVAL(INDEX)*EP

6620 LO=INT(DMIN/INTERVAL)

6630 HI=-INT(-DMAX/INTERVAL)

6640 WHILE (HI-LO)>MAX.INTERVAL

6650　 INDEX=INDEX+1

6660　 IF INDEX>10 THEN INDEX=1:EP=EP*10

6670　 INTERVAL=INTERVAL(INDEX)*EP

6680　 LO=INT(DMIN/INTERVAL)

6690　 HI=-INT(-DMAX/INTERVAL)

6700 WEND

6710 SMIN=LO*INTERVAL

6720 SMAX=HI*INTERVAL

6730 DIF=HI-LO

6740 STP=1

6750 ERASE INTERVAL

6760 RETURN

6770 '

6780 ' ** ラベリング関数 **

6790 '

6800 *SCALE.RTN:

6810 'RE=INT(LOG(ARGUMENT)/LOG(10))

6820 RE=INT(LOG(CDBL(ARGUMENT))/LOG(10#))

6830 'EP=10^RE

6840 EP=VAL("1E"+STR$(RE)):'[EXPONENT PART]

6850 MP=ARGUMENT/EP　　　 :'[MANTISSA PART]

6860 FMP=INT(MP)　:'[FLOORING MANTISSA PART]

6870 CMP=-INT(-MP):'[CEILING MANTISSA PART]

6880 RETURN

6890 '

6900 ' ** 参照間隔 **

6910 '

6920 *REF.INTERVAL:

6930 DATA .1,.2,.5,.5,.5,1,1,1,1,1

6940 '

6950 ' *** REPORT ***

6960 '

6970 *REPORT:

6980 CLS 3

6990 PRINT:PRINT

7000 PRINT "DATA No.=";N1

7010 PRINT

7020 PRINT " "," MEAN"," S.D."," MIN"," MAX"

7030 FOR I=1 TO NV

7040 PRINT "X("+MID$(STR$(I),2)+")",AVE(I),DEV(I),XMIN(I),XMAX(I)

7050 NEXT I

7060 PRINT:PRINT

7070 '

7080 PRINT "[ DEFINED EQUATION ]"

7090 PRINT "Y= ";FORM$

7100 '

7110 GOSUB *CEF.CONFIDENCE

7120 GOSUB *TEMP.STOP

7130 '

7140 CLS 3

7150 PRINT:PRINT "残差"

7160 PRINT " No."," OBSERVED Y"," EXPECTED Y"," RESIDUAL"

7170 FOR I=1 TO N1

7180　 FOR V=1 TO NV:X(V)=XX1(V,I):NEXT V

7190　 GOSUB *DEFINE.FORMULA:YH=FNY(X)

7200　 PRINT I,Y1(I),YH,Y1(I)-YH

7210　 IF (I MOD 20)=0 THEN PRINT

7220 NEXT I

7230 PRINT

7240 GOSUB *TEMP.STOP

7250 RETURN

7260 '

7270 *CEF.CONFIDENCE:

7280 PRINT

7290 PRINT "回帰係数とその９５％信頼区間"

7300 DF1=1:DF2=N1-M

7310 GOSUB *CEF.CONFIDENCE.SUB

7320 PRINT:PRINT

7330 RETURN

7340 '

7350 *CEF.CONFIDENCE.SUB:

7360 FOR J=1 TO M

7370　 P$="P("+MID$(STR$(J),2)+")= "

7380　 PRINT P$;:PRINT P(J);

7390　 IF SD(J)=0 THEN PRINT:GOTO 7420

7400　 GOSUB *CEF.ERROR

7410　 PRINT TAB(30);"S.E.=";SE(J)　　　　　　　　　　　　　　　　　　　　　　　　　　　　　;TAB(50);P(J)-SE(J)*UUT;" - ";P(J)+SE(J)*UUT

7420 NEXT J

7430 RETURN

7440 '

7450 *TEMP.STOP:

7460 PRINT "何かキーを押してください"

7470 WHILE INKEY$="":WEND

7480 RETURN

7490 '

7500 ' *** ERROR of COEFFICIENT ***

7510 '

7520 *CEF.ERROR:

7530 PP=.05 :GOSUB *F.PERCENT

7540 UUT=SQR(UUF)

7550 RETURN

7560 '

7570 ' ** NORMAL DISTRIBUTION **

7580 '

7590 *NORMAL.PERCENT:

7600 XXX=-LOG(4*PP*(1-PP))

7610 UUN=SQR(XXX*(2.06118-5.72622/(XXX+11.6406)))

7620 IF PP>.5 THEN UUN=-UUN

7630 RETURN

7640 '

7650 ' ** T-DISTRIBUTION **

7660 '

7670 *T.PERCENT:

7680 GOSUB *NORMAL.PERCENT

7690 UUN=ABS(UUN)

7700 UU2=UUN*UUN

7710 AA1=(UU2+1)/DF/4

7720 AA2=((5*UU2+16)*UU2+3)/96/DF/DF

7730 AA3=(((3*UU2+19)*UU2+17)*UU2-15)/384/DF/DF/DF

7740 UUT=UUN*(1+AA1+AA2+AA3)

7750 IF PP>.5 THEN UUT=-UUT

7760 RETURN

7770 '

7780 ' ** F-DISTRIBUTION **

7790 '

7800 *F.PERCENT:

7810 IF DF1=1 THEN DF=DF2:PP=PP/2:GOSUB *T.PERCENT　　　　　　　　　　　　　　　　　　　　　　　　:UUF=UUT^2:PP=PP*2:RETURN

7820 IF DF2=1 THEN DF=DF2:PP=(1-PP)/2:GOSUB *T.PERCENT　　　　　　　　　　　　　　　　　　　　　　:UUF=1/UUT^2:PP=1-PP*2:RETURN

7830 IF DF1=2 THEN UUF=DF2/2*(PP^(-2/DF2)-1):RETURN

7840 IF DF2=2 THEN UUF=2/DF1/((1-PP)^(-2/DF1)-1):RETURN

7850 AA=2/9/DF1:BB=2/9/DF2:CC=1-AA:DD=1-BB

7860 GOSUB *NORMAL.PERCENT

7870 UUF=CC*DD+UUN*SQR(CC*CC*BB+DD*DD*AA-AA*BB*UUN*UUN)

7880 UUF=(UUF/(DD*DD-BB*UUN*UUN))^3

7890 RETURN

7900 '

7910 ' ** CHI-SQUARE DISTRIBUTION **

7920 '

7930 *CHI2.PERCENT:

7940 IF DF=1 THEN PP=PP/2:GOSUB *NORMAL.PERCENT　　　　　　　　　　　　　　　　　　　　　　　　　:UUX=UUN*UUN:PP=PP*2:RETURN

7950 IF DF=2 THEN UUX=-2*LOG(PP):RETURN

7960 GOSUB *NORMAL.PERCENT

7970 UUX=2/(9*DF)

7980 UUX=1-UUX+UUN*SQR(UUX)

7990 UUX=UUX^3*DF

8000 RETURN

　

＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝

　

［３］係数変換式と誤差伝播の公式

　

　ここでは，コラム：楕円近似について（その１）で，宿題として遺しておいた

　　ａ1，ａ2，ａ3，ａ4，ａ5　←→　ａ，ｂ，θ，ｘ0，ｙ0

双方向の変換式を紹介します．

　

　自分で式を立ててみれば，この問題の面倒さが実感できますが，

　　ａ1，ａ2，ａ3，ａ4，ａ5　←　　ａ，ｂ，θ，ｘ0，ｙ0

は易しいものの，逆方向の変換

　　ａ1，ａ2，ａ3，ａ4，ａ5　　→　ａ，ｂ，θ，ｘ0，ｙ0

は結構面倒です．以下に係数変換プログラムを掲げますが，9140行〜9240行で，

　　ａ1，ａ2，ａ3，ａ4，ａ5　　→　ａ，ｂ，θ，ｘ0，ｙ0

の係数変換を行っています．

　

9000 '

9010 ' ** 係数変換 **

9020 '

9030 *P2A:

9040 C=COS(P(3)):S=SIN(P(3))

9050 K=1-((P(4)*C+P(5)*S)/P(1))^2-((P(4)*S-P(5)*C)/P(2))^2

9060 A1=(C/P(1))^2+(S/P(2))^2

9070 A2=(1/P(1)^2-1/P(2)^2)*C*S*2

9080 A3=(S/P(1))^2+(C/P(2))^2

9090 A4=-(P(4)*C+P(5)*S)*C*2/P(1)^2-(P(4)*S-P(5)*C)*S*2/P(2)^2

9100 A5=-(P(4)*C+P(5)*S)*S*2/P(1)^2+(P(4)*S-P(5)*C)*C*2/P(2)^2

9110 A(1)=A1/K:A(2)=A2/K:A(3)=A3/K:A(4)=A(4)/K:A(5)=A5/K

9120 RETURN

9130 '

9140 *A2P:

9150 A1=A(1):A2=A(2):A3=A(3):A4=A(4):A5=A(5)

9160 DIF=A1-A3

9170 IF DIF=0 THEN P(3)=PI/4 ELSE P(3)=1/2*ATN(A2/(A1-A3))

9180 C=COS(P(3)):S=SIN(P(3))

9190 K=1/(1+(C*C-S*S)/4*((A4*C+A5*S)^2/(A1*C*C-A3*S*S)-(A4*S-A5*C)^2/(A1*S*S-A3*C*C)))

9200 P(1)=SQR((C*C-S*S)/(A1*C*C-A3*S*S)/K)

9210 P(2)=SQR(-(C*C-S*S)/(A1*S*S-A3*C*C)/K)

9220 P(4)=-(P(1)*P(1)*C*(A4*C+A5*S)+P(2)*P(2)*S*(A4*S-A5*C))/2*K

9230 P(5)=-(P(1)*P(1)*S*(A4*C+A5*S)-P(2)*P(2)*C*(A4*S-A5*C))/2*K

9240 RETURN

　

　この変換式が利用できると，まず最初に

　　ｆ（ｘ，ｙ）＝ａ1ｘ^2 ＋ａ2ｘｙ＋ａ3ｙ^2＋ａ4ｘ＋ａ5ｙ−１

に対して，多変量線形最小２乗法によるあてはめを行い，そのあとで，

　　｛（ｘ−ｘ0）ｃｏｓθ＋（ｙ−ｙ0）ｓｉｎθ｝^2／ａ^2＋｛（ｘ−ｘ0）ｓｉｎθ−（ｙ−ｙ0）ｃｏｓθ｝^2／ｂ^2＝１

の形に直してから描画を行うという２ステップの計算方法も考えられ，実際，その方があてはめ計算が短時間で完了します．

　

　また，線形最小２乗法より高速にあてはめ計算する方法は他には考えられませんから，当初，多変量線形最小２乗法による楕円あてはめを採用しておりました．しかし，その後，楕円の基本情報であるａ，ｂ，θ，ｘ0，ｙ0の推定誤差

　　Δａ，Δｂ，Δθ，Δｘ0，Δｙ0

の計算が要求されたことによって，事情は一変しました．

　

　上に掲げた変換式は

　　ａ ＝ｇ1（ａ1，ａ2，ａ3，ａ4，ａ5）

　　ｂ ＝ｇ2（ａ1，ａ2，ａ3，ａ4，ａ5）

　　θ ＝ｇ3（ａ1，ａ2，ａ3，ａ4，ａ5）

　　ｘ0＝ｇ4（ａ1，ａ2，ａ3，ａ4，ａ5）

　　ｙ0＝ｇ5（ａ1，ａ2，ａ3，ａ4，ａ5）

と書くことができるのですが，誤差伝播の公式を使えば，ａの誤差（Δａ）^2は，

　　（Δａ）^2＝ΣΣ（∂ｇ1／∂ａi）（∂ｇ1／∂ａj）ｃｏｖ（ａi，ａj）

として得られます．Δｂ，Δθ，Δｘ0，Δｙ0についても同様です．

　

　この式でｃｏｖ（ａi，ａj）は多変量線形最小２乗法で計算した分散共分散行列であって，それは既に求められているわけですから，式を見る限り，割合簡単に一方の誤差から他方の誤差が求められそうです．ところが，見るとやるとでは大違い．この微分計算はとても面倒で，間違いなしに計算できる自信はなく，プログラムを作ってみる気になれませんでした．結局，このような理由から，多変量非線形最小２乗法を採用することにしました．

　

　非線形最小２乗法を用いると，あてはめ計算に時間がかかりますが，その代わり，パラメータが楕円の基本情報（半径や中心など）そのものになっていますから，描画が短時間で済みますし，また，パラメータの誤差値を簡単に求めることができるという利点もあるのです．

　

＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝

　

【補】代数曲線のはなし

　

　２変数ｘ，ｙの多項式ｆ（ｘ，ｙ）＝０で定義される曲線を平面代数曲線と呼びます．ｆ（ｘ，ｙ）＝０が２次式の場合，その一般式は，

　　ａｘ^2＋ｈｘｙ＋ｂｙ^2＋ｃｘ＋ｄｙ＋ｅ＝０

のごとく，項数６の多項式として書くことができます．２次曲線には楕円，放物線，双曲線があり，それらは円錐（必ずしも直円錐でなくてよい）を平面で切断したときの切り口として現れる一群の曲線，すなわち円錐曲線です．

　

　同様に，３次曲線とはｆ（ｘ，ｙ）＝０が２変数ｘ，ｙの３次あるいは３次以下の方程式で与えられた曲線です．３次曲線の例としては，ディオクレスのシッソイド（ｘ^3＋ｘｙ^2＝ｙ^2）があげられますが，これは古代ギリシアにおいて立方体倍積問題に用いられた曲線です．また，

　　ｙ＝ｘ^3＋ｘ^2＋ｘ＋１

　　ｙ^3＝ｘｙ^2−２ｘ^2ｙ＋ｙ−３

なども３次曲線で，一般式の項数は１０になります．平面内ｎ次曲線ｆ（ｘ，ｙ）＝０の一般式の項数は，

　　3Ｈn＝n+2Ｃn＝（ｎ＋２）（ｎ＋１）／２

で計算されます．

　

　２次曲線の分類については，３種類の円錐曲線，すなわち楕円，双曲線，放物線になることは既に述べたとおりですが，同じことをもっと高次の曲線・曲面に対して考えるのは自然なことでしょう．３次曲線の分類には，２次曲線とは異なった種類の難解さが要求されましたが，ニュートンはあらゆる場合を考察して，最終的に３次曲線は全部で７８種類が必要であることを示すに至り，さらに３次曲線の一般式が５個の標準形に帰することを示しました．

　

　ニュートンの３次曲線の分類に引き続いて，オイラーは４次平面曲線の分類を企てましたが，可能な場合の数が非常に多いという理由で断念しています．この問題に対する答えは長い間知られていなかったのですが，プリュッカーが１９世紀に４次曲線の１５２の型を数え上げることによって解かれました．

　

　また，一直線上にない３点を通る２次曲線，４点を通る３次曲線はただひとつ存在しますが，それは座標軸の方向が定まっている場合：

　　ｙ＝ａｘ^2＋ｂｘ＋ｃ，ｙ＝ａｘ^3＋ｂｘ^2＋ｃｘ＋ｄ

のようにｙ＝ｆ（ｘ）の場合であって，一般には，平面上の任意の位置にある５点が唯一の円錐曲線を決定します．ニュートンは「プリンキピア」のなかで５点を通る円錐曲線の作図法などを案出しながら壮大な天体力学を展開しています．

　

　ｎ次平面代数曲線の方程式ｆ（ｘ，ｙ）＝０は，（ｎ＋１）（ｎ＋２）／２個の係数をもっていますが，ｆに定数を掛けても曲線は変わりませんから，ｎ次曲線は

　　（ｎ＋１）（ｎ＋２）／２−１＝ｎ（ｎ＋３）／２

個のパラメータに依っていることになります．

　

　そこで，平面内に与えられたｎ（ｎ＋３）／２個の点（ｘi，ｙi）を通るという条件によって曲線を決定するという問題が自然に提起されます．ニュートンはこうした研究を応用して，２次曲線上の５点，３次曲線上の９点が与えられた場合にこれを作図する方法を見いだしたのです．

　

＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝

　

【補】射影変換と代数多様体

　

　放物線，楕円，双曲線はまとめて円錐曲線とも呼ばれますが，２次式で定義されるので，２次曲線ともいいます．そして，無限遠点を導入して，考えている曲線を射影曲線として捉えると，２次曲線はひとつのものとして統一的に考えられるようになります（射影幾何）．なぜなら，違いは無限遠直線の選び方（無限遠直線と交わらない，接する，交わる）にあるだけであって，どれも同種の曲線と考えることができるからです．これらが複素代数多様体として，同型であることをみていくことにしましょう．

　

　楕円ｘ^2／ａ^2＋ｙ^2／ｂ^2＝１は，Ｘ＝ｘ／ａ，Ｙ＝ｙ／ｂという変数変換をすれば，円Ｘ^2＋Ｙ^2＝１になりますから，楕円と円とは代数多様体として同型です．

　

　双曲線ｘ^2／ａ^2−ｙ^2／ｂ^2＝１は，実数Ｒ^2の世界では円や楕円と同型ではありませんが，複素数Ｃ^2の世界で考えれば，Ｘ＝ｘ／ａ，Ｙ＝ｉｙ／ｂｉという変数変換によって，円Ｘ^2＋Ｙ^2＝１になりますから，双曲線と円とは複素代数多様体として同型です．

　

　放物線ｙ＝ｘ^2は，複素代数多様体としても円とは同型になりませんが，射影化という操作，すなわち，放物線に１個の無限遠点を付け加えることによって，射影平面Ｐ^2内の曲線ＹＺ＝Ｘ^2となり，これに対してＹ＝Ｘ'＋ｉＹ'，Ｚ＝Ｘ'−ｉＹ'，Ｘ＝Ｚ'と変数変換すれば，射影化された円Ｘ'^2＋Ｙ'^2＝Ｚ'^2（これは円ｘ^2＋ｙ^2＝１に２個の無限遠点を付け加えたもの）となるので，射影化された放物線と円とは同型です．

　

　また，放物線と直線は１：１対応しますから同型です．したがって，このように複素化し射影化すれば，直線，円，楕円，双曲線，放物線はすべて同型になってしまいます．

　

＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝

　

　しかし，ｘ^3＋ｙ^3＝１のような３次曲線はこれらと同型にはなりません．また，３次曲線は，射影変換を用いれば次のいずれかに変換されます．

　　(1)ｙ^2＝ｘ^3

　　(2)ｙ^2＝ｘ^2（ｘ−１）

　　(3)ｙ^2＝ｘ（ｘ−１）（ｘ−λ）

　

　(1)は「く」の字型曲線で原点で尖点をもちます．(2)は「の」の字型曲線で原点を通ったところでループを描いて自分自身と交差しますから，原点が２重点となります．(3)はループと弓形曲線の２つに分離します．すなわち，(1)(2)は特異点をもち，(3)は非特異です．したがって，滑らかな非特異３次曲線は(3)の標準形に表せます．

　

　特異点を有する(1)(2)は

　　ｙ^2＝ｘ^3　→　（ｔ^3，ｔ^2）

　　ｙ^2＝ｘ^2（ｘ−１）　→　（ｔ^2＋１，ｔ（ｔ^2＋１））

より，曲線上のすべての有理点をパラメトライズすることができます．すなわち有理曲線ですが，それに対して，(3)のように，３次曲線が異なる３根をもつ有理係数の多項式の場合は，楕円曲線と呼ばれる非有理曲線で，２次曲線とは本質的に異なってきます．

　

　これらは特異点による分類といってもよいのですが，射影変換によって互いに写り合う３次曲線は同型とみなされます．そこで，３次曲線のｊ-不変量が定義されます．

　

　非特異３次曲線の標準型：

　　ｙ^2＝ｘ（ｘ−１）（ｘ−λ）

のｊ-不変量は

　　ｊ＝２^8（λ^2−λ＋１）^3／λ^2（λ−１）^2

によって定義されます．λ＝−１のときｊ＝１７２８，λ＝−ζ6（１の６乗根）のときｊ＝０となります．

　

　ｊｰ不変量はモジュラー不変量とも呼ばれ，

　　ｊ（λ）＝ｊ（１−λ）＝ｊ（１／λ）

　＝ｊ（１−１／λ）＝ｊ（１／（１−λ））＝ｊ（λ／（１−λ））

ですから，４個の点｛０，１，λ，∞｝の入れ替えに依存しないinvariantで，最も単純で重要な保型関数と考えられます．

　

　複比を

　　λ＝｛（λ0−λ2）／（λ1−λ2）｝／｛（λ0−λ3）／（λ1−λ3）｝

によって定義すると，λiの順序を変えるとλの値は変わります．すなわち，｛λ0，λ1，λ2，λ3｝からつくられる複比の値は，

　　λ，１−λ，１／（１−λ），１／λ，λ／（λ−１），（λ−１）／λ

の６つのどれかに移ります．

　

　この順序による曖昧さを消すために，λの６つの分数変換の不変式をとって，

　　ｊ＝２^8（λ^2−λ＋１）^3／λ^2（λ−１）^2

とおくのです．複比は一次分数変換で不変であり，ｊもまた射影変換で不変です．（直線上の４点の複比は射影によって不変である）

　

　なお，

　　ｊ（λ）＝ｊ（１−λ）＝ｊ（１／λ）

が成り立てば，あとの等式はこの２つから導かれますから，有理関数

　　（λ^2−λ＋１）^3／λ^2（λ−１）^2

が本質的であって，係数２^8には本質的な意味はありません．実際，

　　（ｘ^2−ｘ＋１）^3／ｘ^2（ｘ−１）^2＝（λ^2−λ＋１）^3／λ^2（λ−１）^2

と，変数ｘの方程式を考えると，

λ^2（λ−１）^2（ｘ^2−ｘ＋１）^3−（λ^2−λ＋１）^3ｘ^2（ｘ−１）^2＝０

はλ≠０，１より，６次方程式となり，

　　λ，１−λ，１／（１−λ），１／λ，λ／（λ−１），（λ−１）／λ

のどれを代入しても成り立ちます．重複が生ずるのは

　　λ^2−λ＋１＝０，λ＝１／２，λ＝−１，λ＝２

の場合に限ります．

　

＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝

　

　ｙ＝ａｘ^3＋ｂｘ^2＋ｃｘ＋ｄという方程式で定まる曲線はおなじみの３次曲線ですが，ｙのところがｙ^2に変わるとワイエルシュトラスの楕円曲線：

　　ｙ^2＝ａｘ^3＋ｂｘ^2＋ｃｘ＋ｄ

になります．ただし，ａ，ｂ，ｃ，ｄは有理数で，右辺の３次式は重根をもたないものと仮定します．楕円曲線をワイエルシュトラス形式に制限しても一般性を失いません．実際，どのような楕円曲線もワイエルシュトラス形式の楕円曲線に双有理的に同値だからです．

　

　また，ｘ^2の項の係数はｘ’＝ｘ＋ｂ／３ａと変数変換（カルダノ変換）することによって簡単に消すことができますから，

　　ｙ^2＝ｘ^3＋ａｘ＋ｂ　　　（４ａ^3＋２７ｂ^2≠０）

を楕円曲線と定義しても構いません．４ａ^3＋２７ｂ^2≠０は重根をもたないための条件です（判別式：Δ＝−（４ａ^3＋２７ｂ^2））．

　

　ワイエルシュトラスの標準形：

　　ｙ^2＝ｘ^3＋ａｘ＋ｂ　　　（２^2ａ^3＋３^3ｂ^2≠０）

のｊ-不変量を計算すると，

　　ｊ＝２^8・３^3ｂ^2／（２^2ａ^3＋３^3ｂ^2）

となります．ｊｰ不変量は，２つの楕円曲線が同じｊｰ不変量をもつかどうかなど，３次曲線を分類する（見分ける）ための指標になっているのです．

　

　４次曲線（項数１５）とか５次（項数２１）以上の高次曲線に対しても射影変換を考えることができます．特異点をもつ３次曲線は適当に座標変換（射影変換）すると(1),(2)のどちらかになりましたが，４次曲線では２０タイプあります．その後，５次曲線は２３０余りのタイプに分類されることが示されましたが，ｎ≧６では複雑すぎてよくわからないようです．

　

＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝

　

【補】位相幾何学的射影平面

　

　ゴム膜でできた長方形の４辺（または２辺）を適当に貼り合わせることを考えてみましょう．そして，長方形を時計回りに順に一周するようなラベルをつけることにします．

　

　すると，１組の対辺だけを向きを保ったまま張り合わせる操作はａ０ａ^(-1)０と書くことができます．ａ０ａ^(-1)０によって円環ができあがります．ａ０ａ０すなわち１組の対辺を向きを逆にして張り合わせる操作によって，メビウスの帯ができあがります．円環の境界は２本の円周ですが，メビウスの帯の境界は１本の円周です．

　

　同様に，ａｂｂ^(-1)ａ^(-1)からは球面，ａｂａ^(-1)ｂ^(-1)からは輪環面（トーラス），ａｂａｂ^(-1)からはクラインの壷，ａｂａｂからは射影平面が得られます．

　

　クラインの壷と射影平面は３次元ユークリッド空間の中では描けません．また，イメージするのは難しいのですが，射影平面はメビウスの帯と円板とを縁に沿って貼り合わせたものと同相です．

　

　円環，球面，輪環面は表から裏に行くことはありませんが，メビウスの帯，クラインの壷，射影平面では縁を越えることなしに曲面の裏側にたどりつきます．こういう面を向きづけ不可能な曲面といいます．

　

　これら６種類はすべて２次元多様体の例です．そして，コンパクトな２次元多様体はどんなものでも，円環，メビウスの帯，球面，輪環面，クラインの壷，射影平面の６つのものを適当に連結和したものと同相になることが知られています．これが２次元多様体の分類定理ですが，３次元以上の多様体に関してはこのような分類定理は存在しません．

　

＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝