数値積分

 投稿者:しばっち  投稿日:2019年 8月22日(木)20時30分49秒
  積分区間[a,b]の数値積分を行います。

  /b
  | f(x)dx
  /a

被積分関数 f(x)=1/xとして

数値積分をリーマン和、台形則、シンプソン則等で行います。

  /x
  | 1/tdt=LOG(x) (x>1)
  /1

LET A=1 !'下限
LET B=2 !'上限
LET N=512 !'分割数
PRINT INTEGRAL(A,B,9999)
PRINT INTEGRAL2(A,B,N)
PRINT INTEGRAL3(A,B,N)
PRINT INTEGRAL4(A,B,N)
PRINT INTEGRAL5(A,B,N)
PRINT INTEGRAL6(A,B,N)
PRINT INTEGRAL7(A,B,N)
PRINT INTEGRAL8(A,B,N)
PRINT INTEGRAL9(A,B,N)
PRINT INTEGRAL10(A,B,N)
PRINT INTEGRAL11(A,B,N)
PRINT INTEGRAL12(A,B,N)
PRINT INTEGRAL13(A,B,N)
PRINT INTEGRAL14(A,B,N)
PRINT INTEGRAL15(A,B,N)
PRINT DE(A,B,N)
PRINT LEGENDRE(A,B)
PRINT CLENSHAW(A,B)
PRINT "----------------------"
PRINT LOG(B) !'真値
END

EXTERNAL FUNCTION FUNC(X)
LET FUNC=1/X
END FUNCTION

EXTERNAL FUNCTION INTEGRAL(A,B,N) !'リーマン和
LET H=(B-A)/(N+1)
FOR J=0 TO N
   LET S=S+H*FUNC(A+H*J)
NEXT J
LET INTEGRAL=S
END FUNCTION

EXTERNAL  FUNCTION INTEGRAL2(A,B,N) !'台形則
DIM R(0 TO 1)
LET R(0)=1/2
LET R(1)=1/2
LET H=(B-A)/N
FOR K=0 TO N-1
   FOR J=0 TO 1
      LET S=S+R(J)*H*FUNC(A+H*(K+J))
   NEXT J
NEXT K
LET INTEGRAL2=S
END FUNCTION

EXTERNAL  FUNCTION INTEGRAL3(A,B,N) !'シンプソン則
DIM R(0 TO 2)
LET R(0)=1/3
LET R(1)=4/3
LET R(2)=1/3
LET H=(B-A)/N/2
FOR K=0 TO N-1
   FOR J=0 TO 2
      LET S=S+R(J)*H*FUNC(A+H*(2*K+J))
   NEXT J
NEXT K
LET INTEGRAL3=S
END FUNCTION

EXTERNAL  FUNCTION INTEGRAL4(A,B,N)
DIM R(0 TO 3)
LET R(0)=3/8
LET R(1)=9/8
LET R(2)=9/8
LET R(3)=3/8
LET H=(B-A)/N/3
FOR K=0 TO N-1
   FOR J=0 TO 3
      LET S=S+R(J)*H*FUNC(A+H*(3*K+J))
   NEXT J
NEXT K
LET INTEGRAL4=S
END FUNCTION

EXTERNAL  FUNCTION INTEGRAL5(A,B,N)
DIM R(0 TO 4)
LET R(0)=14/45
LET R(1)=64/45
LET R(2)=8/15
LET R(3)=64/45
LET R(4)=14/45
LET H=(B-A)/N/4
FOR K=0 TO N-1
   FOR J=0 TO 4
      LET S=S+R(J)*H*FUNC(A+H*(4*K+J))
   NEXT J
NEXT K
LET INTEGRAL5=S
END FUNCTION

EXTERNAL  FUNCTION INTEGRAL6(A,B,N)
DIM R(0 TO 5)
LET R(0)=95/288
LET R(1)=125/96
LET R(2)=125/144
LET R(3)=125/144
LET R(4)=125/96
LET R(5)=95/288
LET H=(B-A)/N/5
FOR K=0 TO N-1
   FOR J=0 TO 5
      LET S=S+R(J)*H*FUNC(A+H*(5*K+J))
   NEXT J
NEXT K
LET INTEGRAL6=S
END FUNCTION

EXTERNAL  FUNCTION INTEGRAL7(A,B,N)
DIM R(0 TO 6)
LET R(0)=41/140
LET R(1)=54/35
LET R(2)=27/140
LET R(3)=68/35
LET R(4)=27/140
LET R(5)=54/35
LET R(6)=41/140
LET H=(B-A)/N/6
FOR K=0 TO N-1
   FOR J=0 TO 6
      LET S=S+R(J)*H*FUNC(A+H*(6*K+J))
   NEXT J
NEXT K
LET INTEGRAL7=S
END FUNCTION

EXTERNAL  FUNCTION INTEGRAL8(A,B,N)
DIM R(0 TO 7)
LET R(0)=5257/17280
LET R(1)=25039/17280
LET R(2)=343/640
LET R(3)=20923/17280
LET R(4)=20923/17280
LET R(5)=343/640
LET R(6)=25039/17280
LET R(7)=5257/17280
LET H=(B-A)/N/7
FOR K=0 TO N-1
   FOR J=0 TO 7
      LET S=S+R(J)*H*FUNC(A+H*(7*K+J))
   NEXT J
NEXT K
LET INTEGRAL8=S
END FUNCTION

EXTERNAL  FUNCTION INTEGRAL9(A,B,N)
DIM R(0 TO 8)
LET R(0)=3956/14175
LET R(1)=23552/14175
LET R(2)=-3712/14175
LET R(3)=41984/14175
LET R(4)=-3632/2835
LET R(5)=41984/14175
LET R(6)=-3712/14175
LET R(7)=23552/14175
LET R(8)=3956/14175
LET H=(B-A)/N/8
FOR K=0 TO N-1
   FOR J=0 TO 8
      LET S=S+R(J)*H*FUNC(A+H*(8*K+J))
   NEXT J
NEXT K
LET INTEGRAL9=S
END FUNCTION

EXTERNAL  FUNCTION INTEGRAL10(A,B,N)
DIM R(0 TO 9)
LET R(0)=25713/89600
LET R(1)=141669/89600
LET R(2)=243/2240
LET R(3)=10881/5600
LET R(4)=26001/44800
LET R(5)=26001/44800
LET R(6)=10881/5600
LET R(7)=243/2240
LET R(8)=141669/89600
LET R(9)=25713/89600
LET H=(B-A)/N/9
FOR K=0 TO N-1
   FOR J=0 TO 9
      LET S=S+R(J)*H*FUNC(A+H*(9*K+J))
   NEXT J
NEXT K
LET INTEGRAL10=S
END FUNCTION

EXTERNAL  FUNCTION INTEGRAL11(A,B,N)
DIM R(0 TO 10)
LET R(0)=80335/299376
LET R(1)=132875/74844
LET R(2)=-80875/99792
LET R(3)=28375/6237
LET R(4)=-24125/5544
LET R(5)=89035/12474
LET R(6)=-24125/5544
LET R(7)=28375/6237
LET R(8)=-80875/99792
LET R(9)=132875/74844
LET R(10)=80335/299376
LET H=(B-A)/N/10
FOR K=0 TO N-1
   FOR J=0 TO 10
      LET S=S+R(J)*H*FUNC(A+H*(10*K+J))
   NEXT J
NEXT K
LET INTEGRAL11=S
END FUNCTION

EXTERNAL  FUNCTION INTEGRAL12(A,B,N)
DIM R(0 TO 11)
LET R(0)=4777223/17418240
LET R(1)=49450643/29030400
LET R(2)=-35608243/87091200
LET R(3)=6166523/1935360
LET R(4)=-17591827/14515200
LET R(5)=28404871/14515200
LET R(6)=28404871/14515200
LET R(7)=-17591827/14515200
LET R(8)=6166523/1935360
LET R(9)=-35608243/87091200
LET R(10)=49450643/29030400
LET R(11)=4777223/17418240
LET H=(B-A)/N/11
FOR K=0 TO N-1
   FOR J=0 TO 11
      LET S=S+R(J)*H*FUNC(A+H*(11*K+J))
   NEXT J
NEXT K
LET INTEGRAL12=S
END FUNCTION

EXTERNAL  FUNCTION INTEGRAL13(A,B,N)
DIM R(0 TO 12)
LET R(0)=1364651/5255250
LET R(1)=150048/79625
LET R(2)=-1264644/875875
LET R(3)=3572512/525525
LET R(4)=-3432753/350350
LET R(5)=14586048/875875
LET R(6)=-2090408/125125
LET R(7)=14586048/875875
LET R(8)=-3432753/350350
LET R(9)=3572512/525525
LET R(10)=-1264644/875875
LET R(11)=150048/79625
LET R(12)=1364651/5255250
LET H=(B-A)/N/12
FOR K=0 TO N-1
   FOR J=0 TO 12
      LET S=S+R(J)*H*FUNC(A+H*(12*K+J))
   NEXT J
NEXT K
LET INTEGRAL13=S
END FUNCTION

EXTERNAL  FUNCTION INTEGRAL14(A,B,N)
DIM R(0 TO 13)
LET R(0)=106364763817/402361344000
LET R(1)=731649485593/402361344000
LET R(2)=-22582626859/22353408000
LET R(3)=144926245243/28740096000
LET R(4)=-78862978129/16094453760
LET R(5)=298542743759/44706816000
LET R(6)=-46704658663/33530112000
LET R(7)=-46704658663/33530112000
LET R(8)=298542743759/44706816000
LET R(9)=-78862978129/16094453760
LET R(10)=144926245243/28740096000
LET R(11)=-22582626859/22353408000
LET R(12)=731649485593/402361344000
LET R(13)=106364763817/402361344000
LET H=(B-A)/N/13
FOR K=0 TO N-1
   FOR J=0 TO 13
      LET S=S+R(J)*H*FUNC(A+H*(13*K+J))
   NEXT J
NEXT K
LET INTEGRAL14=S
END FUNCTION

EXTERNAL  FUNCTION INTEGRAL15(A,B,N)
DIM R(0 TO 14)
LET R(0)=631693279/2501928000
LET R(1)=311056753/156370500
LET R(2)=-5395044599/2501928000
LET R(3)=765940609/78185250
LET R(4)=-46375653541/2501928000
LET R(5)=5525678207/156370500
LET R(6)=-39205297537/833976000
LET R(7)=712193069/13030875
LET R(8)=-39205297537/833976000
LET R(9)=5525678207/156370500
LET R(10)=-46375653541/2501928000
LET R(11)=765940609/78185250
LET R(12)=-5395044599/2501928000
LET R(13)=311056753/156370500
LET R(14)=631693279/2501928000
LET H=(B-A)/N/14
FOR K=0 TO N-1
   FOR J=0 TO 14
      LET S=S+R(J)*H*FUNC(A+H*(14*K+J))
   NEXT J
NEXT K
LET INTEGRAL15=S
END FUNCTION

EXTERNAL  FUNCTION DE(A,B,N) !'二重指数関数法(DE法)
LET U=(B-A)/2
LET V=(B+A)/2
LET H=(B-A)/N
FOR K=-4 TO 4 STEP H
   LET X=Q(K)
   LET S=S+H*FUNC(U*X+V)*U*QQ(K)
NEXT K
LET DE=S

FUNCTION Q(X)
LET Q=TANH(PI/2*SINH(X))
END FUNCTION

FUNCTION QQ(X)
LET QQ=PI/2*COSH(X)/COSH(PI/2*SINH(X))^2
END FUNCTION
END FUNCTION

EXTERNAL  FUNCTION LEGENDRE(A,B) !'ガウス・ルジャンドル則(80次)
LET  U=(B+A)/2
LET  V=(B-A)/2
FOR I=1 TO 80
   READ X,W
   LET  S=S+W*FUNC(U+V*X)*V
NEXT I
LET LEGENDRE=S
DATA -.9995538226516306298800805,1.1449500031869415345441719E-03
DATA -.9976498643982376888994942,2.6635335895126816692935358E-03
DATA -.9942275409656882778920635,4.1803131246948952367393042E-03
DATA -.9892913024997555310265032,5.6909224514031986492691071E-03
DATA -.9828485727386290704182880,7.1929047681173127526755709E-03
DATA -.9749091405857277933856452,8.6839452692608584264094522E-03
DATA -.9654850890437992514522732,1.0161766041103064520831850E-02
DATA -.9545907663436349054934815,1.1624114120797826916466770E-02
DATA -.9422427613098726747522660,1.3068761592401339293786826E-02
DATA -.9284598771724457959530460,1.4493508040509076116962075E-02
DATA -.9132631025717576541647337,1.5896183583725688044902909E-02
DATA -.8966755794387706831943241,1.7274652056269306358584207E-02
DATA -.8787225676782138287037733,1.8626814208299031428735414E-02
DATA -.8594314066631110969771921,1.9950610878141998928891929E-02
DATA -.8388314735802552756166230,2.1244026115782006388710737E-02
DATA -.8169541386814634703711250,2.2505090246332461926221590E-02
DATA -.7938327175046054499486393,2.3731882865930101293192525E-02
DATA -.7695024201350413738656161,2.4922535764115491105117847E-02
DATA -.7440002975835972723165405,2.6075235767565117902968744E-02
DATA -.7173651853620998802540683,2.7188227500486380674418707E-02
DATA -.6896376443420276007712076,2.8259816057276862396753198E-02
DATA -.6608598989861198017359671,2.9288369583267847692767586E-02
DATA -.6310757730468719662479284,3.0272321759557980661220010E-02
DATA -.6003306228297517431547463,3.1210174188114701642442867E-02
DATA -.5686712681227097847254858,3.2100498673487773148056490E-02
DATA -.5361459208971319320198573,3.2941939397645401382836181E-02
DATA -.5028041118887849875936728,3.3733214984611522816675163E-02
DATA -.4686966151705444770360784,3.4473120451753928794364227E-02
DATA -.4338753708317560930623867,3.5160529044747593495526592E-02
DATA -.3983934058819692270243796,3.5794393953416054602861589E-02
DATA -.3623047534994873156190433,3.6373749905835978043964991E-02
DATA -.3256643707477019146191129,3.6897714638276008839150997E-02
DATA -.2885280548845118531091393,3.7365490238730490026705377E-02
DATA -.2509523583922721204931588,3.7776364362001397489774976E-02
DATA -.2129945028576661325723885,3.8129711314477638344206792E-02
DATA -.1747122918326468125593390,3.8424993006959423185212436E-02
DATA -.1361640228091438865592411,3.8661759774076463327077110E-02
DATA -.0974083984415845990632785,3.8839651059051968931774183E-02
DATA -.0585044371524206686289933,3.8958395962769531198625525E-02
DATA -.0195113832567939976543512,3.9017813656306654811280439E-02
DATA  .0195113832567939976543512,3.9017813656306654811280439E-02
DATA  .0585044371524206686289933,3.8958395962769531198625525E-02
DATA  .0974083984415845990632785,3.8839651059051968931774183E-02
DATA  .1361640228091438865592411,3.8661759774076463327077110E-02
DATA  .1747122918326468125593390,3.8424993006959423185212436E-02
DATA  .2129945028576661325723885,3.8129711314477638344206792E-02
DATA  .2509523583922721204931588,3.7776364362001397489774976E-02
DATA  .2885280548845118531091393,3.7365490238730490026705377E-02
DATA  .3256643707477019146191129,3.6897714638276008839150997E-02
DATA  .3623047534994873156190433,3.6373749905835978043964991E-02
DATA  .3983934058819692270243796,3.5794393953416054602861589E-02
DATA  .4338753708317560930623867,3.5160529044747593495526592E-02
DATA  .4686966151705444770360784,3.4473120451753928794364227E-02
DATA  .5028041118887849875936728,3.3733214984611522816675163E-02
DATA  .5361459208971319320198573,3.2941939397645401382836181E-02
DATA  .5686712681227097847254858,3.2100498673487773148056490E-02
DATA  .6003306228297517431547463,3.1210174188114701642442867E-02
DATA  .6310757730468719662479284,3.0272321759557980661220010E-02
DATA  .6608598989861198017359671,2.9288369583267847692767586E-02
DATA  .6896376443420276007712076,2.8259816057276862396753198E-02
DATA  .7173651853620998802540683,2.7188227500486380674418707E-02
DATA  .7440002975835972723165405,2.6075235767565117902968744E-02
DATA  .7695024201350413738656161,2.4922535764115491105117847E-02
DATA  .7938327175046054499486393,2.3731882865930101293192525E-02
DATA  .8169541386814634703711250,2.2505090246332461926221590E-02
DATA  .8388314735802552756166230,2.1244026115782006388710737E-02
DATA  .8594314066631110969771921,1.9950610878141998928891929E-02
DATA  .8787225676782138287037733,1.8626814208299031428735414E-02
DATA  .8966755794387706831943241,1.7274652056269306358584207E-02
DATA  .9132631025717576541647337,1.5896183583725688044902909E-02
DATA  .9284598771724457959530460,1.4493508040509076116962075E-02
DATA  .9422427613098726747522660,1.3068761592401339293786826E-02
DATA  .9545907663436349054934815,1.1624114120797826916466770E-02
DATA  .9654850890437992514522732,1.0161766041103064520831850E-02
DATA  .9749091405857277933856452,8.6839452692608584264094522E-03
DATA  .9828485727386290704182880,7.1929047681173127526755709E-03
DATA  .9892913024997555310265032,5.6909224514031986492691071E-03
DATA  .9942275409656882778920635,4.1803131246948952367393042E-03
DATA  .9976498643982376888994942,2.6635335895126816692935358E-03
DATA  .9995538226516306298800805,1.1449500031869415345441719E-03
END FUNCTION

EXTERNAL  FUNCTION CLENSHAW(A,B) !'クレンショ・カーティス則(51次)
LET  U=(B+A)/2
LET  V=(B-A)/2
FOR I=1 TO 51
   READ X,W
   LET  S=S+W*FUNC(U+V*X)*V
NEXT I
LET CLENSHAW=S
DATA   1.00000000000000000e+00 , 4.00160064025610244e-04
DATA   9.98026728428271562e-01 , 3.85156958477016962e-03
DATA   9.92114701314477831e-01 , 7.90787777387773065e-03
DATA   9.82287250728688681e-01 , 1.17571702393763373e-02
DATA   9.68583161128631119e-01 , 1.56353432816211783e-02
DATA   9.51056516295153572e-01 , 1.94096870342574738e-02
DATA   9.29776485888251404e-01 , 2.31345330553685722e-02
DATA   9.04827052466019528e-01 , 2.67490465016775492e-02
DATA   8.76306680043863587e-01 , 3.02721896767790866e-02
DATA   8.44327925502015079e-01 , 3.36647886957266792e-02
DATA   8.09016994374947424e-01 , 3.69334712479578906e-02
DATA   7.70513242775789231e-01 , 4.00489673231533687e-02
DATA   7.28968627421411523e-01 , 4.30127205858580039e-02
DATA   6.84547105928688674e-01 , 4.58012512695708427e-02
DATA   6.37423989748689710e-01 , 4.84138485266917921e-02
DATA   5.87785252292473129e-01 , 5.08310623057877259e-02
DATA   5.35826794978996618e-01 , 5.30515834900600101e-02
DATA   4.81753674101715275e-01 , 5.50591409946138935e-02
DATA   4.25779291565072649e-01 , 5.68527406895258182e-02
DATA   3.68124552684677959e-01 , 5.84188403643660493e-02
DATA   3.09016994374947424e-01 , 5.97573499442651138e-02
DATA   2.48689887164854788e-01 , 6.08571934325762743e-02
DATA   1.87381314585724631e-01 , 6.17195908592974985e-02
DATA   1.25333233564304245e-01 , 6.23357554958117397e-02
DATA   6.27905195293133761e-02 , 6.27085107726588897e-02
DATA   -1.66091007405623902e-130 , 6.28312135806494038e-02
DATA   -6.27905195293133761e-02 , 6.27085107726588897e-02
DATA   -1.25333233564304245e-01 , 6.23357554958117397e-02
DATA   -1.87381314585724631e-01 , 6.17195908592974985e-02
DATA   -2.48689887164854788e-01 , 6.08571934325762743e-02
DATA   -3.09016994374947424e-01 , 5.97573499442651138e-02
DATA   -3.68124552684677959e-01 , 5.84188403643660493e-02
DATA   -4.25779291565072649e-01 , 5.68527406895258182e-02
DATA   -4.81753674101715275e-01 , 5.50591409946138935e-02
DATA   -5.35826794978996618e-01 , 5.30515834900600101e-02
DATA   -5.87785252292473129e-01 , 5.08310623057877259e-02
DATA   -6.37423989748689710e-01 , 4.84138485266917921e-02
DATA   -6.84547105928688674e-01 , 4.58012512695708427e-02
DATA   -7.28968627421411523e-01 , 4.30127205858580039e-02
DATA   -7.70513242775789231e-01 , 4.00489673231533687e-02
DATA   -8.09016994374947424e-01 , 3.69334712479578906e-02
DATA   -8.44327925502015079e-01 , 3.36647886957266792e-02
DATA   -8.76306680043863587e-01 , 3.02721896767790866e-02
DATA   -9.04827052466019528e-01 , 2.67490465016775492e-02
DATA   -9.29776485888251404e-01 , 2.31345330553685722e-02
DATA   -9.51056516295153572e-01 , 1.94096870342574738e-02
DATA   -9.68583161128631119e-01 , 1.56353432816211783e-02
DATA   -9.82287250728688681e-01 , 1.17571702393763373e-02
DATA   -9.92114701314477831e-01 , 7.90787777387773065e-03
DATA   -9.98026728428271562e-01 , 3.85156958477016962e-03
DATA   -1.00000000000000000e+00 , 4.00160064025610244e-04
END FUNCTION
 

戻る