|
|
> No.1802[元記事へ]
エス・テーさんへのお返事です。
> 今回、十進ベーシックを使って、円：x^2+y^2+Ax+By+C=0 と 放物線：y=px^2+qx+r の交点を
> 求めるプログラムをつくろうとしたのですが、どのようなプログラムをつくればよいのか全くわかりません。
円: x^2+y^2+Ax+By+C=0 ←式1
放物線: y=Px^2+Qx+R ←式2
とすると、
式1に式2を代入して、yを消去する。
x^2+(Px^2+Qx+R)^2+Ax+B(Px^2+Qx+R)+C=0
x^2+((Px^2+Qx)+R)^2+Ax+B(Px^2+Qx+R)+C=0
x^2+(Px^2+Qx)^2+2R(Px^2+Qx)+R^2+Ax+B(Px^2+Qx+R)+C=0
x^2+P^2x^4+2PQx^3+Q^2x^2+2PRx^2+2QRx+R^2+Ax+BPx^2+BQx+BR+C=0
整理して、
P^2x^4+2PQx^3+(Q^2+2PR+BP+1)x^2+(2QR+A+BQ)x+(R^2+BR+C)=0
この4次方程式の実数解は、（実数解でなくてもよい）
代数
・オイラー（Euler）の方法
数値計算
・DKA法
などの解法があります。詳しい説明は、数値計算の文献を参照してください。
ここでは、
数値計算
因数分解して、次数を下げていく
・ニュートン法＋組立除法
で求めてみました。
また、計算結果の確認のため、作図関連のプログラムも組み込んでみました。（半分程度）
!図形と方程式
SET WINDOW -8,8,-8,8 !表示領域を設定する ※調整が必要である
DRAW grid !XY座標
PUBLIC NUMERIC gcCOLOR,gcSTYLE,gcLINESTYLE !描画色、線種
LET gcCOLOR=1 !黒色
LET gcLINESTYLE=1 !実線
!------------------------------ ここまでがサブルーチン
LET A=2 !円x^2+y^2+Ax+By+C=0
LET B=-3
LET C=-20
CALL gcDRAWCIRCLE(A,B,C) !円を描く
LET P=1 !放物線y=Px^2+Qx+R
LET Q=2
LET R=-4
CALL gcDRAWFNC2(P,Q,R, -10,10) !放物線を描く
!ニュートン法による１変数の代数方程式の根
LET N=4 !ｎ次
DIM D(0 TO N) !x^nの係数 d(n)*x^n+d(n-1)*x^(n-1)+ … +d(1)*x+d(0)
LET D(4)=P^2
LET D(3)=2*P*Q
LET D(2)=Q^2+2*P*R+B*P+1
LET D(1)=2*Q*R+A+B*Q
LET D(0)=R^2+B*R+C
!ニュートン法で零点（根）を求める
LET cEps=1e-10 !精度 ※要調整
LET M=N
FOR i=1 TO N !ｎ個
LET Xk=COS((4*i-1)*PI/(4*M+2)) !近似根
LET iter=200 !反復回数
FOR k=1 TO iter
CALL FxdFx(M,D,Xk, f,df) !f(Xk),f'(Xk)の算出
LET WXk=Xk-f/df !Xk+1=Xk-f(Xk)/f'(Xk)
IF ABS(WXk-Xk)<=ABS(Xk)*cEps THEN EXIT FOR !収束すれば終了
LET Xk=WXk !次へ
NEXT k
IF k>iter THEN
PRINT iter;"回では収束しません。"
STOP
END IF
PRINT USING "## ###.###############": i,Xk
CALL gcDRAWPOINT(Xk,P*Xk^2+Q*Xk+R) !交点を描く
DIM QQ(0 TO N) !組立除法で、f(x)=(x-α)Q(x)と因数分解する
CALL poly_divByLin(M,D,Xk, QQ,RR)
MAT D=QQ !copy it
LET M=M-1 !次数を下げる
NEXT i
END
EXTERNAL SUB FxdFx(N,A(),x0, f,df) !関数値f(x0)と微分係数f'(x0)を求める ※n≧1
LET f=A(N)
LET df=f
FOR j=N-1 TO 1 STEP -1 !ホーナー法（組立除法）
LET f=f*x0+A(j)
LET df=df*x0+f
NEXT j
LET f=f*x0+A(0)
END SUB
EXTERNAL SUB poly_divByLin(N,A(),v, Q(),R) !多項式a(x)をx-αで割ったときの商q(x)と余りrを求める
MAT Q=ZER
LET s=0
FOR i=N TO 0 STEP -1 !ホーナーの方法
LET Q(i)=s !※「係数にｘを掛けて次の係数を加える」が組立除法になる
LET s=s*v+A(i) !(…(((An*X+An-1)*X+An-2)*X+An-3)*X+…+A1)*X+A0
NEXT i
LET R=s
END SUB
!作図ルーチン
EXTERNAL SUB gcDRAWPOINT(x,y) !点(x,y)を描く
SET AREA COLOR gcCOLOR
DRAW disk WITH SCALE(0.15)*SHIFT(x,y) !※拡大率0.1は調整が必要である
END SUB
EXTERNAL SUB gcDRAWCIRCLE(A,B,C) !円x^2+y^2+Ax+By+C=0を描く
LET RR=(A^2+B^2)/4-C !判別式
IF RR>=0 THEN
SET LINE COLOR gcCOLOR
SET LINE STYLE gcLINESTYLE
LET CX=-A/2 !中心
LET CY=-B/2
LET R=SQR(RR) !半径
FOR i=0 TO 360 !(x-CX)^2+(y-CY)^2=R^2として描く
PLOT LINES: R*COS(RAD(i))+CX,R*SIN(RAD(i))+CY;
NEXT i
PLOT LINES
ELSE
PRINT "半径が負なので、円が成立しません。"; A;B;C
END IF
END SUB
EXTERNAL SUB gcDRAWFNC2(A,B,C,d1,d2) !2次関数y=Ax^2+Bx+C、x=[d1,d2]を描く
IF A=0 THEN
PRINT "A=0なので、2次関数ではありません。"; A;B;C
ELSE
ASK WINDOW x1,x2,y1,y2
LET x1=MAX(x1,d1)
LET x2=MIN(x2,d2)
SET LINE COLOR gcCOLOR
FOR x=x1 TO x2 STEP 1/2^8 !※折れ線による
PLOT LINES: x,A*x^2+B*x+C;
NEXT x
PLOT LINES
END IF
END SUB
|
|
投稿日：2012年 3月 6日(火)11時49分9秒