解曲線

 投稿者：しばっち  投稿日：2020年 3月29日(日)19時40分30秒
  関数 f(x,y)を導関数としてルンゲクッタ法で解き、その解曲線群を描画します。


LET XS=-5
LET XE=5
LET YS=-5
LET YE=5
SET WINDOW XS,XE,YS,YE
DRAW GRID
LET N=100 !'分割数
LET H=(XE-XS)/N
FOR YY=YS TO XE
   FOR XX=XS TO XE
      FOR I=1 TO 2
         LET Y=YY !'初期値 Y
         LET X=XX !'初期値 X
         WHEN EXCEPTION IN
            FOR J=1 TO N
               LET K1=F(X,Y)
               LET K2=F(X+H/2,Y+H/2*K1)
               LET K3=F(X+H/2,Y+H/2*K2)
               LET K4=F(X+H,Y+H*K3)
               LET Y0=Y+H*(K1+2*K2+2*K3+K4)/6
               LET X0=X+H
               PLOT LINES:X,Y;X0,Y0
               LET X=X0
               LET Y=Y0
            NEXT J
         USE
            PLOT LINES
         END WHEN
         LET H=-H !'符号を逆にして反対側(初期値より前側)を描画する
      NEXT I
   NEXT   XX
NEXT YY
END

EXTERNAL  FUNCTION F(X,Y) !'dy/dx=f(x,y)
!'LET F=Y
!'LET F=Y*Y*X
!'LET F=Y*X
!'LET F=Y*X*X
!'LET F=1/(X*X+1)
!'LET F=Y*Y
!'LET F=Y*Y*Y
LET F=2*X
!'LET F=3*X*X
!'LET F=-Y/(X+1)^2
!'LET F=-X/Y !うまく描けない
END FUNCTION
 		 

解曲線

 投稿者：しばっち  投稿日：2020年 3月29日(日)19時44分8秒
  dx/dt=f1(t,x,y) dy/dt=f2(t,x,y)として2元連立常微分方程式を
ルンゲクッタ法で解き、その解曲線群を描画します。
(※とりあえず同心円が描かれます dy/dx=-X/Y)

LET XS=-5
LET XE=5
LET YS=-5
LET YE=5
SET WINDOW XS,XE,YS,YE
DRAW GRID
LET N=50 !'分割数
LET H=(XE-XS)/N
FOR YY=YS TO YE
   FOR XX=XS TO XE
   ! FOR I=1 TO 2
      LET T=0  !'初期値 T ??? いくつでも変わらない?
      LET X=XX !'初期値 X
      LET Y=YY !'初期値 Y
      WHEN EXCEPTION IN
         FOR J=1 TO N
            LET K1=F1(T,X,Y)
            LET L1=F2(T,X,Y)

            LET K2=F1(T+H/2,X+H/2*K1,Y+H/2*L1)
            LET L2=F2(T+H/2,X+H/2*K1,Y+H/2*L1)

            LET K3=F1(T+H/2,X+H/2*K2,Y+H/2*L2)
            LET L3=F2(T+H/2,X+H/2*K2,Y+H/2*L2)

            LET K4=F1(T+H,X+H*K3,Y+H*L3)
            LET L4=F2(T+H,X+H*K3,Y+H*L3)

            LET T0=T+H
            LET X0=X+H*(K1+2*K2+2*K3+K4)/6
            LET Y0=Y+H*(L1+2*L2+2*L3+L4)/6
            PLOT LINES:X,Y;X0,Y0
            LET X=X0
            LET Y=Y0
            LET T=T0
         NEXT J
      USE
         PLOT LINES
      END WHEN
      ! LET H=-H
      ! NEXT I
   NEXT XX
NEXT YY
END

EXTERNAL  FUNCTION F1(T,X,Y) !'dx/dt=f1(t,x,y)
LET F1=Y
END FUNCTION

EXTERNAL  FUNCTION F2(T,X,Y) !'dy/dt=f2(t,x,y)
LET F2=-X
END FUNCTION
------------------------------------------------------------------------------------------------
LET XS=-5
LET XE=5
LET YS=-5
LET YE=5
SET WINDOW XS,XE,YS,YE
DRAW GRID
LET N=100 !'分割数
LET H=(XE-XS)/N
FOR YY=YS TO YE STEP 3
   FOR XX=XS TO XE STEP 3
      FOR I=1 TO 2
         LET T=0  !'初期値 T
         LET X=XX !'初期値 X
         LET Y=YY !'初期値 Y
         WHEN EXCEPTION IN
            FOR J=1 TO N
               LET K1=F1(T,X,Y)
               LET L1=F2(T,X,Y)

               LET K2=F1(T+H/2,X+H/2*K1,Y+H/2*L1)
               LET L2=F2(T+H/2,X+H/2*K1,Y+H/2*L1)

               LET K3=F1(T+H/2,X+H/2*K2,Y+H/2*L2)
               LET L3=F2(T+H/2,X+H/2*K2,Y+H/2*L2)

               LET K4=F1(T+H,X+H*K3,Y+H*L3)
               LET L4=F2(T+H,X+H*K3,Y+H*L3)

               LET X0=X+H*(K1+2*K2+2*K3+K4)/6
               LET Y0=Y+H*(L1+2*L2+2*L3+L4)/6
               LET T0=T+H

               SET LINE COLOR "RED"
               PLOT LINES:T,X;T0,X0

               SET LINE COLOR "BLUE"
               PLOT LINES:T,Y;T0,Y0

               SET LINE COLOR "GREEN"
               PLOT LINES:X,Y;X0,Y0
               LET X=X0
               LET Y=Y0
               LET T=T0
            NEXT J
         USE
            PLOT LINES
         END WHEN
         LET H=-H
      NEXT  I
   NEXT XX
NEXT YY
END

EXTERNAL  FUNCTION F1(T,X,Y) !'dx/dt=f1(t,x,y)
LET F1=Y
END FUNCTION

EXTERNAL  FUNCTION F2(T,X,Y) !'dy/dt=f2(t,x,y)
LET F2=-1
END FUNCTION
 		 

解曲線

 投稿者：しばっち  投稿日：2020年 3月29日(日)19時45分51秒
  dx/dt=f1(t,x,y,z) dy/dt=f2(t,x,y,z) dz/dt=f3(t,x,y,z)として3元連立常微分方程式を
ルンゲクッタ法で解き、その解曲線群を描画します。

RANDOMIZE
LET XS=-4
LET XE=4
LET YS=-4
LET YE=4
SET WINDOW XS,XE,YS,YE
DRAW GRID
LET N=100 !'分割数
LET H=(XE-XS)/N
FOR L=1 TO 15
   LET XX=INT(RND*10)-5
   LET YY=INT(RND*10)-5
   LET ZZ=INT(RND*10)-5
   FOR I=1 TO 2
      LET X=XX !'初期値 X
      LET Y=YY !'初期値 Y
      LET Z=ZZ !'初期値 Z
      LET T=0
      FOR J=1 TO N
         LET K1=F1(T,X,Y,Z)
         LET L1=F2(T,X,Y,Z)
         LET M1=F3(T,X,Y,Z)

         LET K2=F1(T+H/2,X+H/2*K1,Y+H/2*L1,Z+H/2*M1)
         LET L2=F2(T+H/2,X+H/2*K1,Y+H/2*L1,Z+H/2*M1)
         LET M2=F3(T+H/2,X+H/2*K1,Y+H/2*L1,Z+H/2*M1)

         LET K3=F1(T+H/2,X+H/2*K2,Y+H/2*L2,Z+H/2*M2)
         LET L3=F2(T+H/2,X+H/2*K2,Y+H/2*L2,Z+H/2*M2)
         LET M3=F3(T+H/2,X+H/2*K2,Y+H/2*L2,Z+H/2*M2)

         LET K4=F1(T+H,X+H*K3,Y+H*L3,Z+H*M3)
         LET L4=F2(T+H,X+H*K3,Y+H*L3,Z+H*M3)
         LET M4=F3(T+H,X+H*K3,Y+H*L3,Z+H*M3)

         LET T0=T+H
         LET X0=X+H*(K1+2*K2+2*K3+K4)/6
         LET Y0=Y+H*(L1+2*L2+2*L3+L4)/6
         LET Z0=Z+H*(M1+2*M2+2*M3+M4)/6
         SET LINE COLOR "RED"
         PLOT LINES:T,X;T0,X0
         SET LINE COLOR "BLUE"
         PLOT LINES:T,Y;T0,Y0
         SET LINE COLOR "GREEN"
         PLOT LINES:T,Z;T0,Z0
         LET X=X0
         LET Y=Y0
         LET Z=Z0
         LET T=T0
      NEXT J
      LET H=-H
   NEXT I
NEXT L

FUNCTION F1(T,X,Y,Z) !'dx/dt=f1(t,x,y,z)
   LET F1=X+Y+Z
END FUNCTION

FUNCTION F2(T,X,Y,Z) !'dy/dt=f2(t,x,y,z)
   LET F2=-4*X-3*Y-7*Z
END FUNCTION

FUNCTION F3(T,X,Y,Z) !'dz/dt=f3(t,x,y,z)
   LET F3=2*X+Y+5*Z
END FUNCTION
END
----------------------------------------------------------------------------------------
ローレンツの方程式を3元連立常微分方程式としてルンゲクッタ法で解き
Lorenzアトラクタを3D描画します。

https://ja.wikipedia.org/wiki/ローレンツ方程式


LET XS=-40
LET XE=40
LET YS=-40
LET YE=40
LET ZTH=0          ! z軸のまわりの回転角
LET XTH=0          ! x軸のまわりの回転角初期値
LET YTH=0          ! y軸のまわりの回転角初期値
DIM M(4,4),POINT(4),ROTX(4,4),ROTY(4,4)
SET WINDOW XS,XE,YS,YE
MAT M=ROTATE(ZTH)
MAT ROTX=IDN ! x軸のまわりの回転
LET ROTX(2,2)=COS(XTH)
LET ROTX(2,3)=SIN(XTH)
LET ROTX(3,2)=-SIN(XTH)
LET ROTX(3,3)=COS(XTH)
MAT ROTY=IDN ! y軸のまわりの回転
LET ROTY(1,1)=COS(YTH)
LET ROTY(1,3)=-SIN(YTH)
LET ROTY(3,1)=SIN(YTH)
LET ROTY(3,3)=COS(YTH)
MAT M=M*ROTX*ROTY
LET N=4000
LET H=1/128
LET A=10
LET B=28
LET C=8/3
LET X=1 !'初期値 X
LET Y=1 !'初期値 Y
LET Z=1 !'初期値 Z
LET T=0
FOR J=1 TO N
   LET K1=F1(T,X,Y,Z)
   LET L1=F2(T,X,Y,Z)
   LET M1=F3(T,X,Y,Z)

   LET K2=F1(T+H/2,X+H/2*K1,Y+H/2*L1,Z+H/2*M1)
   LET L2=F2(T+H/2,X+H/2*K1,Y+H/2*L1,Z+H/2*M1)
   LET M2=F3(T+H/2,X+H/2*K1,Y+H/2*L1,Z+H/2*M1)

   LET K3=F1(T+H/2,X+H/2*K2,Y+H/2*L2,Z+H/2*M2)
   LET L3=F2(T+H/2,X+H/2*K2,Y+H/2*L2,Z+H/2*M2)
   LET M3=F3(T+H/2,X+H/2*K2,Y+H/2*L2,Z+H/2*M2)

   LET K4=F1(T+H,X+H*K3,Y+H*L3,Z+H*M3)
   LET L4=F2(T+H,X+H*K3,Y+H*L3,Z+H*M3)
   LET M4=F3(T+H,X+H*K3,Y+H*L3,Z+H*M3)

   LET T=T+H
   LET X=X+H*(K1+2*K2+2*K3+K4)/6
   LET Y=Y+H*(L1+2*L2+2*L3+L4)/6
   LET Z=Z+H*(M1+2*M2+2*M3+M4)/6
   CALL PLOT(X,Y,Z)
NEXT J

SUB PLOT(X,Y,Z)
   LET POINT(1)=X
   LET POINT(2)=Y
   LET POINT(3)=Z
   MAT POINT=POINT*M
   PLOT LINES:POINT(1),POINT(2);
END SUB

FUNCTION F1(T,X,Y,Z) !'dx/dt=f1(t,x,y,z)
   LET F1=A*(Y-X)
END FUNCTION

FUNCTION F2(T,X,Y,Z) !'dy/dt=f2(t,x,y,z)
   LET F2=X*(B-Z)-Y
END FUNCTION

FUNCTION F3(T,X,Y,Z) !'dz/dt=f3(t,x,y,z)
   LET F3=X*Y-C*Z
END FUNCTION
END
 		 

戻る