|
|
!'/1
!'| f(x)dx/(1-x)^α/(1+x)^β
!'/-1
LET N=10
LET ALPHA=1 !'α,β>-1 整数のみ
LET BETA=2
LET H=1/N/8
DIM X(N),W(N)
PRINT " 零点 重み"
FOR XX=-1+H TO 1-H STEP H
LET LO=XX
LET HI=XX+H
IF P(ALPHA,BETA,N,LO)*P(ALPHA,BETA,N,HI)<0 THEN
DO
LET M=(LO+HI)/2
IF P(ALPHA,BETA,N,LO)*P(ALPHA,BETA,N,M)<0 THEN LET HI=M ELSE LET LO=M
LOOP UNTIL ABS(HI-LO)<1E-13
LET I=I+1
LET X(I)=(LO+HI)/2
LET W(I)=LAMDA(ALPHA,BETA,N,X(I))
PRINT I;":";X(I);W(I)
END IF
NEXT XX
INPUT PROMPT "下限 =":A
INPUT PROMPT "上限 =":B
LET U=(B+A)/2
LET V=(B-A)/2
FOR I=1 TO N
LET S=S+W(I)*F(U+V*X(I))*V/(1-X(I))^ALPHA/(1+X(I))^BETA
NEXT I
PRINT S
PRINT B^3/3-A^3/3
END
EXTERNAL FUNCTION F(X)
LET F=X*X
END FUNCTION
EXTERNAL FUNCTION P(A,B,N,Z) !'Jacobi polynomials ※ガンマ関数gamma(n+1)を階乗fact(n)で代用
FOR M=0 TO N
LET S=S+COMB(N,M)*FACT(A+B+N+M)/FACT(A+M)*((Z-1)/2)^M
NEXT M
LET P=S*FACT(A+N)/FACT(N)/FACT(A+B+N)
END FUNCTION
EXTERNAL FUNCTION DIFFP(K,A,B,N,Z) !'d^k/dz^k P(A,B,N,Z)
LET DIFFP=FACT(A+B+N+K)/2^K/FACT(A+B+N)*P(A+K,B+K,N-K,Z)
END FUNCTION
EXTERNAL FUNCTION LAMDA(A,B,N,X)
LET LAMDA=-(2*N+A+B+2)/(N+A+B+1)*FACT(N+A)*FACT(N+B)/FACT(N+A+B)/FACT(N+1)*2^(A+B)/DIFFP(1,A,B,N,X)/P(A,B,N+1,X)
END FUNCTION
|
|