|
投稿者：しばっち
投稿日：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
|
|
|