オイラー法

 投稿者:しばっち  投稿日:2020年 2月16日(日)12時09分39秒
  https://ja.wikipedia.org/wiki/オイラー法
https://ja.wikipedia.org/wiki/ルンゲ=クッタ法
http://www.yamamo10.jp/yamamoto/lecture/2005/5E/test_2/html/node3.html
http://www.caero.mech.tohoku.ac.jp/publicData/Daiguji/Chapter7.pdf
https://www.ktech.biz/jp/article/10-ode-1/
https://www.ktech.biz/jp/article/10-ode-2/

!'  y'=F(x,y)
OPTION BASE 0
LET XS=0 !' XS~XEまで
LET XE=1
LET X=XS
LET Y=0 !'初期値
LET N=65536 !'分割数
LET H=(XE-XS)/N
LET NN=8
DIM XX(NN),YY(NN),A(NN)
FOR I=1 TO N
   LET Y=Y+F(X,Y)*H
   LET X=X+H
   IF MOD(I,N/NN)=0 THEN
      LET K=K+1
      LET XX(K)=X
      LET YY(K)=Y
      PRINT X;Y;X^2
   END IF
NEXT I
CALL CALC(NN,XX,YY,A)
CALL DISPLAY(NN-1,A)
END

EXTERNAL  FUNCTION F(X,Y) !' y'=f(x,y)
LET F=2*X
END FUNCTION

EXTERNAL  SUB CALC(N,XX(),YY(),A())
OPTION BASE 0
DIM C(N),X(N),Y(N)
REDIM A(N)
LET X(1)=1
FOR I=1 TO N
   MAT Y=ZER
   LET Y(0)=1
   FOR L=1 TO N
      MAT C=ZER
      IF I<>L THEN
         LET X(0)=-XX(L)
         FOR J=0 TO N
            FOR K=0 TO N
               IF J+K<=N THEN LET C(K+J)=C(K+J)+X(K)*Y(J)
            NEXT  K
         NEXT  J
         MAT Y=C
      END IF
   NEXT  L
   LET H=1
   FOR J=1 TO N
      IF I<>J THEN
         LET H=H*(XX(I)-XX(J))
      END IF
   NEXT J
   LET H=YY(I)/H
   MAT Y=(H)*Y
   MAT A=A+Y
NEXT I
END SUB

EXTERNAL SUB DISPLAY(N,A())
FOR I=0 TO N
   IF ABS(A(I))<1E-5 THEN
      LET A(I)=0
   ELSE
      LET A(I)=NUM(A(I))
   END IF
NEXT I
IF N>1 THEN
   IF A(I)<>0 THEN
      LET FLG=1
      IF A(N)<0 THEN PRINT "-";
      IF ABS(A(N))<>1 THEN
         PRINT STR$(ABS(A(N)));"*X^";STR$(N);
      ELSE
         PRINT "X^";STR$(N);
      END IF
   END IF
END IF
FOR I=N-1 TO 2 STEP-1
   IF A(I)<>0 THEN
      IF A(I)<0 THEN
         PRINT "-";
         LET FLG=1
      ELSE
         IF FLG=1 THEN PRINT "+";
      END IF
      IF ABS(A(I))<>1 THEN
         LET FLG=1
         PRINT STR$(ABS(A(I)));"*X^";STR$(I);
      ELSEIF ABS(A(I))=1 THEN
         LET FLG=1
         PRINT "X^";STR$(I);
      END IF
   END IF
NEXT I
IF A(1)<>0 THEN
   IF N>1 THEN
      IF A(1)<0 THEN
         PRINT "-";
         LET FLG=1
      ELSE
         IF FLG=1 THEN PRINT "+";
      END IF
   END IF
   IF ABS(A(1))<>1 THEN
      LET FLG=1
      PRINT STR$(ABS(A(1)));"*X";
   ELSEIF ABS(A(1))=1 THEN
      LET FLG=1
      PRINT "X";
   END IF
END IF
IF A(0)<>0 THEN
   IF A(0)<0 THEN
      PRINT "-";
   ELSE
      IF FLG=1 THEN PRINT "+";
   END IF
   PRINT STR$(ABS(A(0)));
END IF
PRINT
END SUB

EXTERNAL FUNCTION NUM(X)
LET EPS=1E-5
FOR I=0 TO 4
   FOR J=1 TO 99
      FOR K=0 TO 1
         LET XX=J*ABS(X)*10^I+K*EPS
         IF ABS(XX)-INT(ABS(XX))<EPS THEN
            LET NUM=SGN(X)*INT(ABS(XX))/10^I/J
            EXIT FUNCTION
         END IF
      NEXT K
   NEXT J
NEXT I
LET NUM=X
END FUNCTION
 

戻る