teacup. [ 掲示板 ] [ 掲示板作成 ] [ 有料掲示板 ] [ ブログ ]

 投稿者
  題名
  内容 入力補助画像・ファイル<IMG>タグが利用可能です。(詳細)
    
 URL
[ ケータイで使う ] [ BBSティッカー ] [ 書込み通知 ] [ 検索 ]

スレッド一覧

  1. スレッドが使えます(2)
  2. Paract BASIC(20)
  3. Amusement_Program(10)
  4. 改修予定のJIS非互換(3)
  5. 複数ページ長編プログラム(新規投稿)(16)
  6. 「十進BASIC第2掲示板」投稿記事リスト(17)
  7. Full BASIC互換ライブラリ(8)
  8. 「十進BASIC掲示板過去ログ」インデックス(トピック)(17)
  9. 人の色覚の数理(14)
  10. 「十進BASIC掲示板過去ログ」インデックス(ツリー)(91)
スレッド一覧(全10)  他のスレッドを探す  スレッド作成

*掲示板をお持ちでない方へ、まずは掲示板を作成しましょう。無料掲示板作成


[018]周期的境界条件下の電子状態(2次元)

 投稿者:mike  投稿日:2017年 6月26日(月)09時33分32秒
返信・引用
   2次元-周期的境界条件下の下で定常状態の電子のエネルギーと電子状態を求めるプログラム
(018harmonicsPSD2D.bas)を公開します。これは、最急降下法による2次元の定常状態の電子状態を求めるプログラム
[011]に周期的境界条件を付けたものです。

表示の説明:
 2次元-周期的境界条件下の自由空間(V(x,y)=0)では基底状態は|0>=constant, E(0)=0 になります。
x、y方向の周期が同じ場合、E1=E2=E3=E4で縮退しているため、|1>..|4>はこの固有空間の中の
任意の直交する関数の組でよいため、形が一意には決まりません。
E5=E6=E7=E8も縮退しています。

試験環境:
本プログラムは十進BASIC 6.6.3.1 / macOS 10.7.5, 十進BASIC Ver 7.8.0 / windows 10 でテストしました。


------------

!
! ========= periodic steepest descent method 2D ==========
!
! 018periodicPSD2D.bas
!   Copyright(C) 2017 Mitsuru Ikeuchi
!   Released under the MIT license ( https://opensource.org/licenses/MIT )
!
! ver 0.0.1   2017.06.26  created
!
OPTION ARITHMETIC NATIVE
DECLARE EXTERNAL SUB psd2d.setInitialCondition,psd2d.SDiteration,psd2d.drawState
DECLARE EXTERNAL FUNCTION INKEY$

LET stateMax = 10  !state 0,1,...,stateMax-1
LET vIndex = 0     !0:free space,  1:2d-crystal
LET iterMax = 5    !5 -iteration in SDiteration()
LET drawMode = 1   !0:draw3DPsi, 1:drawPsi
CALL setInitialCondition(stateMax,vIndex)

DO
   CALL SDiteration(stateMax, iterMax)
   CALL drawState(drawMode,vIndex)
   LET S$=INKEY$
   IF S$="6" THEN CALL setInitialCondition(stateMax,0) !Hermonics
   IF S$="7" THEN CALL setInitialCondition(stateMax,1) !Well
   IF S$="8" THEN LET drawMode = 0 !0:draw3DPsi
   IF S$="9" THEN LET drawMode = 1 !1:drawPsi
LOOP UNTIL S$=CHR$(27)
END

EXTERNAL FUNCTION INKEY$  !from decimal BASIC library
OPTION ARITHMETIC NATIVE
SET ECHO "OFF"
LET S$=""
CHARACTER INPUT NOWAIT: S$
LET INKEY$=S$
END FUNCTION

! ---------- periodic steepest descent method 2D ----------
!
! system Hamiltonian: H = -delta/2 + V(r) , delta r = div grad r
! eigen energy Ei, eigen function set { |i> }
!
! procedure : successive approximation
! (i) trial function set { |0>,|1>,..,|i>,.. }
! (2) energy of |i> : ei = <i|H|i>/<i|i>
! (3) steepest gradient direction (H-ei)|i>
! (4) next generation : |i(next)> = |i> - dampingFactor*(H-ei)|i>
! (5) orthogonalization { |0>,|1>,..,|i>,.. }  (Gram-Schmidt)
! (6) goto (2)
!
MODULE psd2d
MODULE OPTION ARITHMETIC NATIVE
OPTION BASE 0
PUBLIC SUB setInitialCondition, SDiteration, drawState
SHARE NUMERIC NNx, NNy, dx, dy, iterCount
SHARE NUMERIC cosAx,sinAx,cosAy,sinAy,cx0,cy0,cz0 !--- 3D graphics
SHARE NUMERIC sdEnergy(20)        ! electron state energy
SHARE NUMERIC sdState(20,200,200) ! electron states 0...20 0:ground state
SHARE NUMERIC wrk(200,200)        ! state work space in steepestDescent
SHARE NUMERIC vv(200,200)         ! external potential
SHARE NUMERIC psicol(64,64)       ! MAT PLOT CELL matrix for psi(x,y)
SHARE NUMERIC vcol(64,64)         ! MAT PLOT CELL matrix  for potential V(x,y)
LET NNx = 64                      ! max number of sdState(,NNx,NNy)
LET NNy = 64                      ! max number of sdState(,NNx,NNy)
LET dx = 0.25                     ! (au) x division
LET dy = 0.25                     ! (au) y division
LET iterCount = 0                 ! sd iteration count

! ---------- set initial condition

EXTERNAL SUB setInitialCondition(stateMax,vIndex)   !public
   DECLARE EXTERNAL SUB setInitialState,setPotential,initDraw
   LET iterCount = 0
   CALL setInitialState(stateMax)
   CALL setPotential(vIndex)
   CALL initDraw
   ! set window
   LET xMargin = 60
   LET yMargin = 120
   SET WINDOW 0,500, 0,500
END SUB

EXTERNAL SUB setInitialState(stateMax)
   DECLARE EXTERNAL SUB normalizeState
   RANDOMIZE
   FOR ist=0 TO stateMax-1
      FOR i=0 TO NNx-1
         FOR j=0 TO NNx-1
            LET sdState(ist,i,j) = RND-0.5
         NEXT j
      NEXT i
      CALL normalizeState(ist)
   NEXT ist
END SUB

EXTERNAL SUB setPotential(vIndex)
   LET x0 = 0.5*NNx*dx
   LET y0 = 0.5*NNy*dy
   FOR i=0 TO NNx-1
      FOR j=0 TO NNy-1
         LET x = i*dx
         LET y = j*dy
         IF vIndex=0 THEN
            LET vv(i,j) = 0
         ELSEIF vIndex=1 THEN
            IF (x-x0)*(y-y0)>0 THEN LET vv(i,j)=10 ELSE LET vv(i,j)=0
         ELSE
            LET vv(i,j)=0
         END IF
      NEXT j
   NEXT i
END SUB

EXTERNAL SUB initDraw !--- set set color pallet and MAT PLOT matrix vcol(,)
   FOR i = 0 TO 50 !--- set color pallet
      SET COLOR MIX( 40+i) 0.02*i,0.02*i,0     ! red for |psi> +
      SET COLOR MIX(100+i) 0,0.02*i,0.02*i     ! blue for |psi> -
      SET COLOR MIX(160+i) 0,0.02*i,0     ! green for V(x)
   NEXT i
   FOR i=0 TO NNx-1 !--- set vcol(,)
      FOR j=0 TO NNy-1
         LET col = 0.02*vv(i,j)
         IF col>1 THEN LET col = 1
         IF col<0 THEN LET col = 0
         LET vcol(i,j) = 160+INT(col*50)
      NEXT j
   NEXT i
END SUB

! ---------- steepest descent iteration

EXTERNAL SUB SDiteration(stateMax, iterMax) !public
   DECLARE EXTERNAL FUNCTION steepestDescent
   DECLARE EXTERNAL SUB GramSchmidt,sortState
   LET damp = 0.01 !damping factor in steepest descent
   FOR i=0 TO iterMax-1
      FOR ist=0 TO stateMax-1
         LET sdEnergy(ist) = steepestDescent(ist, damp)
      NEXT ist
      CALL GramSchmidt(stateMax)
      CALL sortState(stateMax)
      LET iterCount = iterCount + 1
   NEXT i
END SUB

EXTERNAL FUNCTION steepestDescent(ist,damp)  !---  steepest descent method
   DECLARE EXTERNAL FUNCTION energyOfState
   DECLARE EXTERNAL SUB normalizeState
   LET h2 = 2*dx*dx
   LET ei = energyOfState(ist)               !---  E_ist = <ist|H|ist>
   FOR i=0 TO NNx-1                          !---  |wrk> = (H-E_ist)|ist>
      LET ip1 = MOD(i+1,NNx)
      LET im1 = MOD(i-1+NNx,NNx)
      FOR j=0 TO NNy-1
         LET jp1 = MOD(j+1,NNy)
         LET jm1 = MOD(j-1+NNy,NNy)
         LET pij = sdState(ist,i,j)
         LET wrk(i,j) = (4*pij-sdState(ist,ip1,j)-sdState(ist,im1,j)-sdState(ist,i,jm1)-sdState(ist,i,jp1))/h2+(vv(i,j)-ei)*pij
      NEXT j
   NEXT i
   FOR i=0 TO NNx-1                          !---  |ist(next)> = |ist> - damp*|wrk>  ( norm(|ist(next)>) <>1 )
      FOR j=0 TO NNy-1
         LET sdState(ist,i,j) = sdState(ist,i,j)-damp*wrk(i,j)
      NEXT j
   NEXT i
   CALL normalizeState(ist)
   LET steepestDescent = ei
END FUNCTION

EXTERNAL FUNCTION energyOfState(ist) !--- E_ist = <ist|H|ist>/<ist|ist>
   LET h2 = 2*dx*dx
   LET s = 0
   LET sn=0
   FOR i=0 TO NNx-1
      LET ip1 = MOD(i+1,NNx)
      LET im1 = MOD(i-1+NNx,NNx)
      FOR j=0 TO NNy-1
         LET jp1 = MOD(j+1,NNy)
         LET jm1 = MOD(j-1+NNy,NNy)
         LET pij = sdState(ist,i,j)
         LET pij = sdState(ist,i,j)
         LET s = s+pij*((4*pij-sdState(ist,ip1,j)-sdState(ist,im1,j)-sdState(ist,i,jm1)-sdState(ist,i,jp1))/h2+vv(i,j)*pij)
         LET sn = sn+pij*pij
      NEXT j
   next i
   LET energyOfState = s/sn
END FUNCTION

EXTERNAL SUB GramSchmidt(stateMax)  !--- Gram-Schmidt orthonormalization
   DECLARE EXTERNAL FUNCTION innerProduct
   DECLARE EXTERNAL SUB normalizeState
   CALL normalizeState(0)
   FOR istate=1 TO stateMax-1
      FOR ist=0 TO istate-1
         LET s = innerProduct(ist,istate)
         FOR i=0 TO NNx-1
            FOR j=0 TO NNy-1
               LET sdState(istate,i,j) = sdState(istate,i,j) - s*sdState(ist,i,j)
            NEXT j
         NEXT i
      NEXT ist
      CALL normalizeState(istate)
   NEXT iState
END SUB

EXTERNAL SUB sortState(stateMax)
   FOR ist=stateMax-2 TO 0 STEP -1
      IF sdEnergy(ist)>sdEnergy(ist+1)+0.00001 THEN
         FOR i=0 TO NNx-1
            FOR j=0 TO NNy-1
               LET w = sdState(ist,i,j)
               LET sdState(ist,i,j) = sdState(ist+1,i,j)
               LET sdState(ist+1,i,j) = w
            NEXT j
         NEXT i
         LET w = sdEnergy(ist)
         LET sdEnergy(ist) = sdEnergy(ist+1)
         LET sdEnergy(ist+1) = w
      END IF
   NEXT ist
END SUB

EXTERNAL FUNCTION innerProduct(ist,jst)  !--- <ist|jst>
   LET s = 0
   FOR i=0 TO NNx-1
      FOR j=0 TO NNy-1
         LET s = s + sdState(ist,i,j)*sdState(jst,i,j)
      NEXT j
   NEXT i
   LET innerProduct = s*dx*dy
END FUNCTION

EXTERNAL SUB normalizeState(ist)
   LET s = 0
   FOR i=0 TO NNx-1
      FOR j=0 TO NNy-1
         LET s = s + sdState(ist,i,j)*sdState(ist,i,j)*dx*dy
      NEXT j
   NEXT i
   LET a = SQR(1/s)
   FOR i=0 TO NNx-1
      FOR j=0 TO NNy-1
         LET sdState(ist,i,j) = a*sdState(ist,i,j)
      NEXT j
   NEXT i
END SUB

! ---------- utility

!EXTERNAL FUNCTION iterationCount
!   LET iterationCount = iterCount
!END FUNCTION

EXTERNAL FUNCTION stateEnergy(ist)
   LET stateEnergy = sdEnergy(ist)
END FUNCTION

! ---------- 3D graphics

EXTERNAL SUB setRotateXYParameters(angleX,angleY,xCenter,yCenter,zCenter)
   LET cosAx = COS(angleX)
   LET sinAx = SIN(angleX)
   LET cosAy = COS(angleY)
   LET sinAy = SIN(angleY)
   LET cx0 = xCenter
   LET cy0 = yCenter
   LET cz0 = zCenter
END SUB

EXTERNAL SUB plotLines3D(x,y,z,xShift,yShift) !shift*xRotateAx*yRotateAy*(shift^-1)
   LET x1 = cosAy*(x-cx0)+sinAy*(sinAx*(y-cy0)+cosAx*(z-cz0))
   LET y1 = cosAx*(y-cy0)-sinAx*(z-cz0)
   LET z1 =-sinAy*(x-cx0)+cosAy*(sinAx*(y-cy0)+cosAx*(z-cz0))
   PLOT LINES: x1+cx0+xShift, y1+cy0+yShift; !z=z1+cz0
END SUB

! ---------- drawState

EXTERNAL SUB drawState(drawMode,vIndex)
   DECLARE EXTERNAL SUB setRotateXYParameters,drawState3D,draw3DPsix6,drawPsix6
   IF vIndex=0 THEN LET v$="free space" ELSE LET v$="2D-crystal"
   SET DRAW MODE HIDDEN
   CLEAR
   IF drawMode=0 THEN
      CALL setRotateXYParameters(-PI/6,-PI/12,NNx/2,NNy/2,0)
      CALL draw3DPsix6
   ELSEIF drawMode=1 THEN
      CALL drawPsix6
   END IF
   !--- caption
   SET TEXT HEIGHT 10
   SET TEXT COLOR 1 ! black
   PLOT TEXT, AT 50,100 ,USING "iteration count =##### ":iterCount
   PLOT TEXT, AT 50, 85 ,USING "E0 =###.###########  E5 =###.###########(au)":sdEnergy(0),sdEnergy(5)
   PLOT TEXT, AT 50, 70 ,USING "E1 =###.###########  E6 =###.###########(au)":sdEnergy(1),sdEnergy(6)
   PLOT TEXT, AT 50, 55 ,USING "E2 =###.###########  E7 =###.###########(au)":sdEnergy(2),sdEnergy(7)
   PLOT TEXT, AT 50, 40 ,USING "E3 =###.###########  E8 =###.###########(au)":sdEnergy(3),sdEnergy(8)
   PLOT TEXT, AT 50, 25 ,USING "E4 =###.###########  E9 =###.###########(au)":sdEnergy(4),sdEnergy(9)
   PLOT TEXT, AT 50, 10 :"periodic steepest descent method 2D"
   PLOT TEXT, AT 50,480 :"potential -> [6] free space [7] 2d-crystal "
   PLOT TEXT, AT 50,465 :"display   -> [8] psi3D      [9] psi  "
   PLOT TEXT, AT 50,450 :"[esc] exit "
   PLOT TEXT, AT 50,420 :"potential : "&v$
   SET DRAW MODE EXPLICIT
END SUB

EXTERNAL SUB draw3DPsix6
   DECLARE EXTERNAL SUB drawState3D
   CALL drawState3D(0,2,0.5, 30,120) !(ist,sc,zMag,xShift,yShift)
   CALL drawState3D(1,2,0.5,180,120)
   CALL drawState3D(2,2,0.5,330,120)
   CALL drawState3D(3,2,0.5, 30,270)
   CALL drawState3D(4,2,0.5,180,270)
   CALL drawState3D(5,2,0.5,330,270)
   SET TEXT HEIGHT 10
   PLOT TEXT, AT  30,130:"|0>"
   PLOT TEXT, AT 180,130:"|1>"
   PLOT TEXT, AT 330,130:"|2>"
   PLOT TEXT, AT  30,280:"|3>"
   PLOT TEXT, AT 180,280:"|4>"
   PLOT TEXT, AT 330,280:"|5>"
END SUB

EXTERNAL SUB drawState3D(ist,sc,zMag,xShift,yShift)
   DECLARE EXTERNAL SUB plotLines3D
   FOR j=0 TO NNy-1 STEP 1
      FOR i=0 TO NNx-1
         LET psi = sdState(ist,i,j)*200
         IF psi>1 THEN
            SET LINE COLOR 4 ! red
         ELSEIF psi<-1 THEN
            SET LINE COLOR 2 ! blue
         ELSE
            SET LINE COLOR 3 ! potential:green
         END IF
         CALL plotLines3D(i*sc,j*sc,zMag*(psi + 2*vv(i,j)),xShift,yShift) !(x,y,z)
      NEXT i
      PLOT LINES
   NEXT j
   FOR i=0 TO NNx-1 STEP 1
      FOR j=0 TO NNy-1
         LET psi = sdState(ist,i,j)*200
         IF psi>1 THEN
            SET LINE COLOR 4 ! red
         ELSEIF psi<-1 THEN
            SET LINE COLOR 2 ! blue
         ELSE
            SET LINE COLOR 3 ! potential:green
         END IF
         CALL plotLines3D(i*sc,j*sc,zMag*(psi + 2*vv(i,j)),xShift,yShift) !(x,y,z)
      NEXT j
      PLOT LINES
   NEXT i
END SUB

EXTERNAL SUB drawPsix6
   DECLARE EXTERNAL SUB drawPsi
   CALL drawPsi(0, 30,135) !(ist,xPos,yPos)
   CALL drawPsi(1,180,135)
   CALL drawPsi(2,330,135)
   CALL drawPsi(3, 30,285)
   CALL drawPsi(4,180,285)
   CALL drawPsi(5,330,285)
   SET TEXT HEIGHT 10
   PLOT TEXT, AT  30,120:"|0>"
   PLOT TEXT, AT 180,120:"|1>"
   PLOT TEXT, AT 330,120:"|2>"
   PLOT TEXT, AT  30,270:"|3>"
   PLOT TEXT, AT 180,270:"|4>"
   PLOT TEXT, AT 330,270:"|5>"
END SUB

EXTERNAL SUB drawPsi(ist,xPos,yPos)  !drawMode=0
   MAT psicol = vcol !--- psicol(,) <- Vcol(,) setted SUB initDraw
   FOR i=0 TO NNx-1 !--- set psicol(,) for MAT PLOT CELLS
      FOR j=0 TO NNy-1
         LET psi = sdState(ist,i,j)
         IF abs(psi)>0.002 THEN
            IF psi>=0 THEN
               LET col = psi*10
               IF col>1 THEN LET col = 1
               IF col<0 THEN LET col = 0
               LET psicol(i,j) = 40+INT(col*50)
            ELSE
               LET col = -psi*10
               IF col>1 THEN LET col = 1
               IF col<0 THEN LET col = 0
               LET psicol(i,j) = 100+INT(col*50)
            END IF
         END IF
      NEXT j
   NEXT i
   MAT PLOT CELLS,IN xPos,yPos; xPos+128,yPos+128 :psicol
END SUB

END MODULE
 
 

Re: 実行時内部エラーの報告

 投稿者:SHIRAISHI Kazuo  投稿日:2017年 6月24日(土)08時49分45秒
返信・引用
  十進BASIC32ビット版は,手続き呼び出し時にCPUのSPレジスタを操作して受け渡し場所を確保しているのですが,それがMAC-OSのお気に召さなかったようです。
 

Re: 実行時内部エラーの報告

 投稿者:SHIRAISHI Kazuo  投稿日:2017年 6月23日(金)23時01分21秒
返信・引用
  ご報告ありがとうございました。

Windows版,Linux版では問題なく,Macだけの問題のようです。

 

実行時内部エラーの報告

 投稿者:mike  投稿日:2017年 6月23日(金)10時41分53秒
返信・引用
  実行時内部エラーが発生しました。

access violation エラーです。

decimalBASIC 6.6.3.1 / MacOS 10.7.5 / MACmini (intel core-i7 2.7GHz 4GB)

再起動後に実行してもエラーが再現します。

下記のプログラムです。
コメントアウトして発生場所を絞り込むと、
SUB setInitialCondition で発生しているものと推定できます。

   !set density
   FOR i=0 TO NNp-1
      LET s = 0.0
      FOR j=0 TO NNp-1
         LET s = s + mass(j)*kernel(distance(i,j), (hh(i)+hh(j))/2.0)
      NEXT j
      LET density(i) = s
   NEXT i

の部分を下記のように変更するとエラーは回避できそうですが

   !set density
   FOR i=0 TO NNp-1
      LET s = 0.0
      FOR j=0 TO NNp-1
         LET r = distance(i,j)
         LET s = s + mass(j)*kernel(r, (hh(i)+hh(j))/2.0)
      NEXT j
      LET density(i) = s
   NEXT i


まだ未完ですが、全体のプログラムを載せます。

-------------

!
! ========= smoothed particle hydrodynamics 2D ==========
!
! 0xxgasSPH2D.bas
!
!    Copyright(C) 2017 Mitsuru Ikeuchi
!    Released under the MIT license ( https://opensource.org/licenses/MIT )
!
!    ver 0.0.0  2017.06.22 created
!
OPTION ARITHMETIC NATIVE
DECLARE EXTERNAL SUB sph2d.setInitialCondition, sph2d.timeEvolution, sph2d.drawParticles

CALL setInitialCondition

LET t0 = TIME
FOR it=1 TO 1000
!CALL timeEvolution(10)
!CALL drawParticles
NEXT it
PRINT TIME - t0
END

! ---------- smoothed particle hydrodynamics 2D ----------
!
! - particle base Lagrangian method
!    Monaghan; "Smoothed Particle Hydrodynamics"
!              Annu. Rev. Astron. Astrophys.1992. 30:543-74
!
!    W(x-xi,h) : kernel weight function (q=|x-xi|/h)
!      = aKernel*(1-1.5*q^2+0.75*q^3)  (q<1)
!      = aKernel*0.25*(2-q)^3  (1<=q<2)
!      = 0  (q>=2)
!
!        aKernel = 2/3/h (1D)
!                = 10/(7Pi)/h^2 (2D)
!                = 1/Pi/h^3 (3D)
!
!    f(x) --> sum[mj/rhoj*f(xj)*W(x-xj,h), j]
!    grad f(x) --> -sum[mj/rhoj*f(xj)*grad W(x-xj,h), j]
!
! - time step
!   (1) registration - set section  (not implemented in this code)
!   (2) set density
!   (3) set pressure
!   (4) Verlet step1 (t = t + dt/2)
!   (5) set acceleration ax ay and power <-- x(t+dt/2),y(t+dt/2)
!   (6) Verlet step2 (t = (t + dt/2) + dt/2)
!   (7) set boundary
!   goto (1)
!
MODULE sph2d
MODULE OPTION ARITHMETIC NATIVE
MODULE OPTION BASE 0
PUBLIC SUB setInitialCondition,timeEvolution,drawParticles
SHARE NUMERIC sysTime, dt, xMax, yMax, hh0, NNp
SHARE NUMERIC xx(5000),yy(5000)  ! (xx(i),yy(i))  : position of i-th particle
SHARE NUMERIC vx(5000),vy(5000)  ! (vx(i),vy(i))  : velocity of i-th particle
SHARE NUMERIC ax(5000),ay(5000)  ! (ax(i),ay(i))  : total force/mass of i-th particle
SHARE NUMERIC hh(5000),mass(5000),cp(5000)
SHARE NUMERIC density(5000),pressure(5000),energy(5000),power(5000)
LET sysTime = 0.0  ! (s) system time
LET dt = 0.0001    ! (s) time step
LET xMax = 5.0     ! (m) x-Box size
LET yMax = 5.0     ! (m) y-Box size
LET hh0 = 0.25     ! (m) h of weight function
LET NNp = 400      ! number of smoothed particle

! ---------- set initial condition

EXTERNAL SUB setInitialCondition
!  i,j, s
   DECLARE EXTERNAL FUNCTION kernel,distance
   RANDOMIZE
   LET sysTime = 0.0
   !set smoothed particles
   FOR i=0 TO NNp-1
      LET mass(i) = (0.029/0.0224)/4.0 !air
      LET cp(i) = 1000.0 !air
      LET xx(i) = 0.1*MOD(i,20)+0.5
      LET yy(i) = 0.1*INT(i/20)+1.0
      LET vx(i) = 0.0
      LET vy(i) = 0.0
      LET ax(i) = 0.0
      LET ay(i) = 0.0
      LET hh(i) = hh0
   NEXT i
   ! set energy and power
   FOR i=0 TO NNp-1
      LET energy(i) = mass(i)*cp(i)*300 ! 300=temp(K)
      LET power(i) = 0.0
   NEXT i
   !set density
   FOR i=0 TO NNp-1
      LET s = 0.0
      FOR j=0 TO NNp-1
         LET s = s + mass(j)*kernel(distance(i,j), (hh(i)+hh(j))/2.0)
      NEXT j
      LET density(i) = s
   NEXT i
   ! set window
   SET WINDOW 0,500, 0,500
END SUB

! ---------- time evolution

EXTERNAL SUB timeEvolution(nCalc)
!  i
   DECLARE EXTERNAL SUB timeStep
   FOR i=0 TO NNp-1
      CALL timeStep
      LET sysTime = sysTime + dt
   NEXT i
END SUB

EXTERNAL SUB timeStep
!  i,j, dtv2,s,rr
   DECLARE EXTERNAL SUB accCalc
   DECLARE EXTERNAL FUNCTION kernel,distance
   LET dtv2 = dt/2.0

   !(2) set density
   FOR i=0 TO NNp-1
      LET s = 0.0
      FOR j=0 TO NNp-1
         LET s = s + mass(j)*kernel(distance(i,j), (hh(i)+hh(j))/2.0)
      Next j
      LET density(i) = s
   NEXT i

   !(3) SET pressure
   FOR i=0 TO NNp-1
      LET pressure(i) = 8.31*density(i)*(energy(i)/(mass(i)*cp(i)))
   NEXT i

   !(4) Verlet step1 (t = t + dt/2)
   FOR i=0 TO NNp-1
      LET vx(i) = vx(i) + dtv2*ax(i)
      LET vy(i) = vy(i) + dtv2*ay(i)
      LET xx(i) = xx(i) + vx(i)*dt
      LET yy(i) = yy(i) + vy(i)*dt
   NEXT i
   !(5) set acceleration ax[] ay[] and power[]
   CALL accCalc
   !(6) Verlet step2 (t = t + dt/2 + dt/2)
   FOR i=0 TO NNp-1
      LET vx(i) = vx(i) + dtv2*ax(i)
      LET vy(i) = vy(i) + dtv2*ay(i)
      LET energy(i) = energy(i) + power(i)*dt
   NEXT i

   !(7) boundary
   LET rr = 1.0
   FOR i=0 TO NNp-1
      if (xx(i) < 0.0) then
         LET xx(i) = 0.0
         LET vx(i) = -rr*vx(i)
         LET vy(i) = rr*vy(i)
      END IF
      if (xx(i) > xMax) then
         LET xx(i) = xMax
         LET vx(i) = -rr*vx(i)
         LET vy(i) = rr*vy(i)
      end if
      if (yy(i) < 0.0) then
         LET yy(i) = 0.0
         LET vx(i) = rr*vx(i)
         LET vy(i) = -rr*vy(i)
      END IF
      if (yy(i) > yMax) then
         LET yy(i) = yMax
         LET vx(i) = rr*vx(i)
         LET vy(i) = -rr*vy(i)
      END IF
   NEXT i
END SUB

! --- (5) set acceleration ax[] ay[] and power[]

EXTERNAL SUB accCalc
!  i,j, ai,aj,b,rij,gradW,gradWx,gradWy
   DECLARE EXTERNAL FUNCTION dwwvdr
   FOR i=0 TO NNp-1
      LET ax(i) = 0.0
      LET ay(i) = 0.0
      LET power(i) = 0.0
      LET ai = pressure(i)/(density(i)*density(i))
      FOR j=0 TO NNp-1
         LET rij = SQR((xx(i)-xx(j))*(xx(i)-xx(j))+(yy(i)-yy(j))*(yy(i)-yy(j)))
         IF (rij>0) THEN
            LET aj = pressure(j)/(density(j)*density(j))
            LET b = mass(j)*(ai+aj)
            LET gradW = dwwvdr(rij, (hh(i)+hh(j))/2.0)
            LET gradWx = gradW*(xx(i)-xx(j))/rij
            LET gradWy = gradW*(yy(i)-yy(j))/rij
            LET ax(i) = ax(i) - b*gradWx
            LET ay(i) = ay(i) - b*gradWy
            LET power(i) = power(i) + 0.5*b*((vx(i)-vx(j))*gradWx+(vy(i)-vy(j))*gradWy)
         END IF
      NEXT j
   NEXT i
END SUB

! --- smoothed particle

EXTERNAL FUNCTION distance(i,j)
!  x,y
   LET x = xx(i)-xx(j)
   LET y = yy(i)-yy(j)
   LET distance = SQR(x*x+y*y)
END FUNCTION

EXTERNAL FUNCTION kernel(r,h)
!  a,q,ret
   LET ret = 0.0
   LET q = r/h
   LET a = 10.0/(7.0*3.141592)/(h*h)
   if (q<1.0) then
      LET ret = a*(1.0-1.5*q*q+0.75*q*q*q)
   ELSEIF (q<2.0) THEN
      LET ret = a*0.25*(2.0-q)*(2.0-q)*(2.0-q)
   END IF
   LET kernel = ret
END FUNCTION

EXTERNAL FUNCTION dwwvdr(r,h)
!  a,q,ret
   LET ret = 0.0
   LET q = r/h
   LET a = 10.0/(7.0*3.141592)/(h*h*h)
   if (q<1.0) then
      LET ret = a*(-3.0*q+2.25*q*q)
   elseif (q<=2.0) then
      LET ret = -a*0.75*(2.0-q)*(2.0-q)
   END IF
   LET dwwvdr = ret
END FUNCTION

! ---------- utility



! ---------- draw particles

EXTERNAL SUB drawParticles
   DECLARE EXTERNAL FUNCTION systemTemprature
   LET boxsize = 300
   LET sc = boxsize/xMax
   LET xp = 100
   LET yp = 80
   LET vScate = 100*dt
   LET fScale = 1000*dt*dt/mass(0)

   SET DRAW MODE HIDDEN
   CLEAR
   SET LINE COLOR 1 ! black : wall
   PLOT LINES: xp,yp; xp+boxSize,yp; boxSize+xp,boxSize+yp; xp,boxSize+yp; xp,yp
   !draw Balls, velocity and force
   FOR i=1 TO nMolec
      SET LINE COLOR 2 ! blue : molecule
      DRAW circle WITH SCALE(sigma/2.0*sc)*SHIFT(xx(i)*sc+xp,yy(i)*sc+yp)
      SET LINE COLOR 4 ! red : velocity
      PLOT LINES: xx(i)*sc+xp,yy(i)*sc+yp;
      PLOT LINES: (xx(i)+vx(i)*vScate)*sc+xp,(yy(i)+vy(i)*vScate)*sc+yp
      !SET LINE COLOR 1 ! black : force
      !PLOT LINES: xx(i)*sc+xp,yy(i)*sc+yp;
      !PLOT LINES: (xx(i)+ffx(i)*fScale)*sc+xp,(yy(i)+ffy(i)*fScale)*sc+yp
   NEXT i
   SET TEXT HEIGHT 6
   SET LINE COLOR 4 ! red : velocity
   PLOT LINES: xp+0.6*boxsize,yp+boxsize+25;xp+0.6*boxsize+50,yp+boxsize+25
   SET TEXT COLOR 4 ! red
   PLOT TEXT, AT xp+0.6*boxsize+55,yp+boxsize+22: "velosity"
   SET LINE COLOR 1 ! black : force
   PLOT LINES: xp+0.6*boxsize,yp+boxsize+10;xp+0.6*boxsize+50,yp+boxsize+10
   SET TEXT COLOR 1 ! black
   PLOT TEXT, AT xp+0.6*boxsize+55,yp+boxsize+7: "force"

   SET TEXT HEIGHT 10
   SET TEXT COLOR 1 ! black
   PLOT TEXT, AT 50, 50 ,USING "time =####.#### ":sysTime
   PLOT TEXT, AT 50, 35 ,USING  "N =####":NNp
   PLOT TEXT, AT 50, 20 :"Ar in the box (2-dimensional molecular dynamics)"
   SET TEXT HEIGHT 6
   PLOT TEXT, AT xp,boxSize+2+yp ,USING  "box size =##.# x ##.# (m)":xMax,yMax
   SET DRAW MODE EXPLICIT
END SUB

END MODULE
 

[017]異種金属結晶の接触

 投稿者:mike  投稿日:2017年 6月19日(月)10時12分40秒
返信・引用
  Morseポテンシャルによる異種金属結晶の接触の分子動力学法プログラム017mixMMD2D.basを公開します。
異なった種類のナノ結晶を接触させます。接触すると表面エネルギーを解放し、温度が高くなります。
一般に表面エネルギー(Morseポテンシャルの場合は解離エネルギーに相当)の大きい方の表面に
小さい方の原子が覆う形が安定になります。

表示の説明:
鉄原子は暗い水色、アルミニウム原子は茶色の円であらわします。
表示は[D]のキーを押すことで変更できます。

試験環境:
本プログラムは十進BASIC 6.6.3.1 / macOS 10.7.5, 十進BASIC Ver 7.8.0 / windows 10 でテストしました。


---------------------


!
! ========= molecular dynamics 2D - Morse potential ==========
!
! 017mixMMD2D.bas
!   Copyright(C) 2017 Mitsuru Ikeuchi
!   Released under the MIT license ( https://opensource.org/licenses/MIT )
!
! ver 0.0.0  2017.04.22  created
! ver 0.0.1  2017.06.11  bug fixed
!
OPTION ARITHMETIC NATIVE
DECLARE EXTERNAL SUB mmd2d.setInitialCondition, mmd2d.moveParticles, mmd2d.drawParticles
DECLARE EXTERNAL FUNCTION INKEY$
LET tempMode = 0    !tempMode: 0:adiabatic  1:constant temperature
LET contTemp = 300  !contTemp: controled constant temperature(K)
LET drawMode = 0    !drawMode: 0:bond  1:circle+bond 2:velocitySpace
DIM menu$(8)
LET menu$(1) = "Fe Al"
LET menu$(2) = "Fe Mo"
LET menu$(3) = "Fe Hg"
LET menu$(4) = "W Pb"
LET menu$(5) = "Ni Cr"
LET menu$(6) = "Ag Cu"
LET menu$(7) = "Al Ar"
LET menu$(8) = "continue"

!setInitialCondition(kind1,kind2,boxSizeInNM,contTemp)
!molKind: 0:W 1:Mo 2:Cr 3:Fe 4:Ni 5:Al 6:Pb 7:Cu 8:Ag 17:Ar 20:Hg
CALL setInitialCondition(3,5,6.0,contTemp)

DO
   CALL moveParticles(tempMode,contTemp)
   CALL drawParticles(tempMode,contTemp,drawMode)
   LET S$=INKEY$
   IF S$="0" THEN
      LOCATE VALUE (1),RANGE 0 TO 4000 : temp
      IF temp>1 THEN
         LET tempMode = 1
         LET contTemp = temp
      ELSE
         LET tempMode = 0
      END IF
   ELSEIF S$="D" OR S$="d" THEN
      LET drawMode = MOD(drawMode+1,3)
   ELSEIF S$="9" THEN
      LOCATE CHOICE (menu$) :nmenu
      IF nmenu=1 THEN      !--- Fe Al box=6x6nm
         CALL setInitialCondition(3,5,6.0,contTemp)
      ELSEIF nmenu=2 THEN  !--- Fe Mo
         CALL setInitialCondition(3,1,6.0,contTemp)
      ELSEIF nmenu=3 THEN  !--- Fe Hg
         CALL setInitialCondition(3,20,6.0,contTemp)
      ELSEIF nmenu=4 THEN  !--- W Pb
         CALL setInitialCondition(0,6,6.0,contTemp)
      ELSEIF nmenu=5 THEN  !--- Ni Cr
         CALL setInitialCondition(4,2,6.0,contTemp)
      ELSEIF nmenu=6 THEN  !--- Ag Cu
         CALL setInitialCondition(8,7,6.0,contTemp)
      ELSEIF nmenu=7 THEN  !--- Al Ar
         CALL setInitialCondition(5,17,6.0,contTemp)
      ELSEIF nmenu=8 THEN  !contine
      END IF
   END IF
LOOP UNTIL S$=CHR$(27)
END


EXTERNAL FUNCTION INKEY$  !from decimalBASIC library
OPTION ARITHMETIC NATIVE
SET ECHO "OFF"
LET S$=""
CHARACTER INPUT NOWAIT: S$
LET INKEY$=S$
END FUNCTION

! ---------- molecular dynamics 2D - Morse potential ----------
!
!  method: velocity Verlet algorithm
!    (1) vi = vi + (Fi/mi)*(0.5dt)
!    (2) ri = ri + vi*dt
!    (3) calculation Fi <- {r1,r2,...,rn} Fi=sum(Fij,j=1 to n),Fij=F(ri-rj)
!    (4) vi = vi + (Fi/mi)*(0.5dt)
!    (6) goto (1)
!  potential: Morse V(r) = D*((1-EXP(-A*(r-r0)))^2-1)
!                        = D*(EXP(-2*A*(r-r0))-2*EXP(-A*(r-r0)))
!    (D:dissociation energy, r0:bond length, A:width parameter { A=SQR(k/(2*D)) }
!             force F(r) = -dV(r)/dr
!                        = 2*D*A*y*(y-1), y=EXP(-A*(r-r0))
MODULE mmd2d
MODULE OPTION ARITHMETIC NATIVE
PUBLIC SUB setInitialCondition !(molKind,boxSizeInNM,xtalSizeInNM,contTemp)
PUBLIC SUB moveParticles !(tempMode,contTemp)
PUBLIC SUB drawParticles !(drawMode:0:realSpace 1:velocitySpace)
SHARE NUMERIC molecKind1,molecKind2, sysTime, dt, nMolec, xMax, yMax, rCutoff, hh
SHARE NUMERIC xx(500),yy(500)           ! (xx(i),yy(i))  : position of i-th particle
SHARE NUMERIC vx(500),vy(500)           ! (vx(i),vy(i))  : velocity of i-th particle
SHARE NUMERIC ffx(500),ffy(500)         ! (ffx(i),ffy(i)): total force of i-th particle
SHARE NUMERIC kind(500),mas(500)        ! kind(i),mas(i) : molec kind, mass of i-th particle
SHARE NUMERIC reg(500,0 TO 100)         ! register near molec reg(i,0):number of near i-th molec
SHARE NUMERIC molecData(0 TO 20,0 TO 3) ! molecule 0:mass, 1:epsilon, 2:sigma, 3:dt
SHARE STRING molecSTR$(0 TO 20)         ! molec string
SHARE NUMERIC forceTable(0 TO 20, 0 TO 20,0 TO 1001) ! force table

LET molecKind1 = 3         ! molecule kind1 - 3:Fe
LET molecKind2 = 1         ! molecule kind2 - 1:Mo
LET sysTime = 0.0          ! system time (s) in the module
LET dt = 5.0*1.0e-15       ! time step (s)
LET nMolec = 100           ! number of particles
LET xMax = 6.0E-9          ! x-Box size (m)
LET yMax = 6.0E-9          ! y-Box size (m)
LET rCutoff = 1.0e-9       ! force cutoff radius (m)
LET hh = 1e-12             ! force table step (m)
LET molecData(0,0)=0       ! if molecData(0,0)=0 then set molecData(,)

! ---------- set initial condition

EXTERNAL SUB setInitialCondition(kind1,kind2,boxSizeInNM,contTemp)
   DECLARE EXTERNAL SUB setMolecData,setMolecules,ajustVelocity,setForceTable
   DECLARE EXTERNAL FUNCTION setCrystalBlock
   RANDOMIZE
   CALL setMolecData
   LET molecKind1 = kind1
   LET molecKind2 = kind2
   LET sysTime = 0.0
   LET xMax = boxSizeInNM*1.0e-9
   LET yMax = boxSizeInNM*1.0e-9
   ! set particles
   LET s = 0.1*xMax
   LET xtalSize = 0.35*xMax
   LET ss = 0.55*xMax
   LET ii = setCrystalBlock(1, kind1, s, s, xtalSize, xtalSize, PI/4)
   LET ii = setCrystalBlock(ii+1, kind2, s, ss, xtalSize, xtalSize, 0)
   LET ii = setCrystalBlock(ii+1, kind2, ss, s, xtalSize, xtalSize, 0)
   LET nMolec = setCrystalBlock(ii+1, kind1, ss, ss, xtalSize, xtalSize, 0)
   CALL ajustVelocity(contTemp)
   LET rCutoff = MIN(1.0e-9, 3.0*MAX(molecData(kind1,3),molecData(kind2,3)))
   CALL setForceTable
   SET COLOR MIX(240) 1.0,0.2,0.2 !for bond direction color
   SET COLOR MIX(241) 0.6,0.6,0.0
   SET COLOR MIX(242) 0.0,1.0,0.0
   SET COLOR MIX(243) 0.0,0.6,0.6
   SET COLOR MIX(244) 0.2,0.2,1.0
   SET COLOR MIX(245) 0.8,0.0,0.8
   ! set window
   SET WINDOW 0,500,0,500
END SUB

EXTERNAL SUB setMolecData
! Morse potential data
!       0:mass(AU) 1:D(eV)  2:A(m^-1)  3:r0(m)
   DATA   183.85 , 0.9906, 1.4116e10, 3.032e-10 !  0 W
   DATA    95.94 , 0.8032, 1.5079e10, 2.976e-10 !  1 Mo
   DATA    51.996, 0.4414, 1.5721e10, 2.754e-10 !  2 Cr
   DATA    55.847, 0.4174, 1.3885e10, 2.845e-10 !  3 Fe
   DATA    58.71 , 0.4205, 1.4199e10, 2.780e-10 !  4 Ni
   DATA    26.98 , 0.2703, 1.1646e10, 3.253e-10 !  5 Al
   DATA   207.19 , 0.2348, 1.1836e10, 3.733e-10 !  6 Pb
   DATA    63.54 , 0.3429, 1.3588e10, 2.866e-10 !  7 Cu
   DATA   107.87 , 0.3323, 1.3690e10, 3.115e-10 !  8 Ag
   DATA    40.08 , 0.1623, 0.8054e10, 4.569e-10 !  9 Ca
   DATA    87.62 , 0.1513, 0.7878e10, 4.988e-10 ! 10 Sr
   DATA   137.34 , 0.1416, 0.6570e10, 5.373e-10 ! 11 Ba
   DATA    22.99 , 0.0633, 0.5900e10, 5.336e-10 ! 12 Na
   DATA    39.102, 0.0542, 0.4977e10, 6.369e-10 ! 13 K
   DATA    85.47 , 0.0464, 0.4298e10, 7.207e-10 ! 14 Rb
   DATA   132.905, 0.0449, 0.4157e10, 7.557e-10 ! 15 Cs
   DATA    20.183, 0.0031, 1.6500e10, 3.076e-10 ! 16 Ne
   DATA    39.948, 0.0104, 1.3400e10, 3.816e-10 ! 17 Ar
   DATA    83.80 , 0.0141, 1.2500e10, 4.097e-10 ! 18 Kr
   DATA   131.30 , 0.0200, 1.2400e10, 4.467e-10 ! 19 Xe
   DATA   200.59 , 0.0734, 1.4900e10, 3.255e-10 ! 20 Hg
   DATA "W" ,"Mo","Cr","Fe","Ni","Al","Pb","Cu","Ag","Ca","Sr"
   DATA "Ba","Na","K" ,"Rb","Cs","Ne","Ar","Kr","Xe","Hg"

   IF molecData(0,0)=0 THEN
      MAT READ molecData
      FOR i=0 TO 20
         LET molecData(i,0) = molecData(i,0)*1.661e-27 !mass(AU)--> (kg)
         LET molecData(i,1) = molecData(i,1)*1.602e-19 !eps(eV) --> (J)
      NEXT i
      MAT READ molecSTR$
   END IF
END SUB

EXTERNAL SUB setGas(nMolecule)
   DECLARE EXTERNAL SUB ajustVelocity
   LET r0 = molecData(knd,3)
   FOR j=1 TO nMolecule
      LET loopCount = 0
      DO
         LET xx(j) = (xMax-2*r0)*RND + r0
         LET yy(j) = (yMax-2*r0)*RND + r0
         FOR i=1 TO j-1
            IF (xx(i)-xx(j))^2+(yy(i)-yy(j))^2 < 2*r0^2 THEN EXIT FOR
         NEXT i
         LET loopCount = loopCount + 1
         IF loopCount>1000 THEN EXIT DO
      LOOP UNTIL i>=j
      IF loopCount>1000 THEN
         LET nMolec = j - 1
         EXIT FOR
      END IF
   NEXT j
   FOR i=1 TO nMolec
      LET vx(i) = 200.0*(RND+RND+RND+RND+RND+RND-3)
      LET vy(i) = 200.0*(RND+RND+RND+RND+RND+RND-3)
      LET ffx(i) = 0.0
      LET ffy(i) = 0.0
   NEXT i
END SUB

EXTERNAL FUNCTION setCrystalBlock(ii, knd, x0, y0, xLen, yLen, theta)
   DECLARE EXTERNAL SUB setParticle
   LET iip = ii
   LET a = 0.98*molecData(knd,3) !r0 of knd
   LET b = 0.866025*a
   LET leng = xLen
   IF (leng<yLen) THEN LET leng = yLen
   LET leng = 1.5*leng
   LET nx = INT(leng/b) + 1
   LET ny = INT(leng/a) + 1
   LET sth = SIN(theta)
   LET cth = COS(theta)
   FOR i=0 TO nx-1
      LET x = b*i - leng/2.0
      FOR j=0 TO ny-1
         LET y = a*j - leng/2.0
         IF MOD(i,2)=1 THEN LET y = y + 0.5*a
         LET xp = x0 + xLen/2.0 + cth*x - sth*y
         LET yp = y0 + yLen/2.0 + sth*x + cth*y
         IF (xp>=x0 AND xp<=x0+xLen AND yp>=y0 AND yp<=y0+yLen) THEN
            CALL setParticle(iip, knd, xp, yp)
            LET iip = iip + 1
         END IF
      NEXT j
   NEXT i
   LET setCrystalBlock = iip - 1
end function

EXTERNAL SUB setParticle(i, knd, x, y)
   LET xx(i) = x
   LET yy(i) = y
   LET vx(i) = 200.0*(RND+RND+RND+RND+RND+RND-3)
   LET vy(i) = 200.0*(RND+RND+RND+RND+RND+RND-3)
   LET ffx(i) = 0.0
   LET ffy(i) = 0.0
   LET mas(i) = molecData(knd,0)
   LET kind(i) = knd
END SUB

EXTERNAL SUB setForceTable
   DECLARE EXTERNAL  FUNCTION cutoff
   FOR ki=0 TO 20
      FOR kj=0 TO 20
         LET dd = SQR(molecData(ki,1)*molecData(kj,1))
         LET aa = 0.5*(molecData(ki,2)+molecData(kj,2))
         LET r0 = 0.5*(molecData(ki,3)+molecData(kj,3))
         FOR ir=10 TO 1001
            LET r = ir*hh
            LET y = EXP(-aa*(r-r0))
            LET forceTable(ki,kj,ir) = cutoff(r)*2.0*dd*aa*y*(y-1)
         NEXT ir
         FOR ir=0 TO 9
            LET forceTable(ki,kj,ir) = forceTable(ki,kj,10)
         NEXT ir
      NEXT kj
   NEXT ki
END SUB

EXTERNAL FUNCTION cutoff(r)
   IF r>0 AND r<0.8*rCutoff THEN
      LET ret = 1
   ELSEIF r>=0.8*rCutoff AND r<rCutoff THEN
      LET ret = 0.5+0.5*COS(PI*(r-0.8*rCutoff)/(0.2*rCutoff))
   else
      LET ret = 0
   END IF
   LET cutoff = ret
END FUNCTION

! ---------- move particles

EXTERNAL SUB moveParticles(tempMode,contTemp) !tempMode 0:adiabatic  1:constant-temp
   DECLARE EXTERNAL SUB moveParticlesDT,ajustVelocity,registerNearMolec
   IF (tempMode=1) THEN CALL ajustVelocity(contTemp)
   CALL registerNearMolec
   FOR i=1 TO 20
      CALL moveParticlesDT
   NEXT i
END SUB

EXTERNAL SUB moveParticlesDT ! velocity Verlet method
   DECLARE EXTERNAL SUB calcForce
   FOR i=1 TO nMolec
      LET a = 0.5*dt/mas(i)
      LET vx(i) = vx(i)+a*ffx(i)
      LET vy(i) = vy(i)+a*ffy(i)
      LET xx(i) = xx(i)+vx(i)*dt
      LET yy(i) = yy(i)+vy(i)*dt
   NEXT i
   CALL calcForce
   FOR i=1 TO nMolec
      LET a = 0.5*dt/mas(i)
      LET vx(i) = vx(i)+a*ffx(i)
      LET vy(i) = vy(i)+a*ffy(i)
   NEXT i
   LET sysTime=sysTime+dt
END SUB

EXTERNAL SUB calcForce
   DECLARE EXTERNAL FUNCTION force,boundaryForce
   LET s = 0.5*3.418e-10
   FOR i=1 TO nMolec
      LET ffx(i) = 0.0
      LET ffy(i) = 0.0
   NEXT i
   FOR i=1 TO nMolec-1
      FOR k=1 TO reg(i,0)-1
         LET j = reg(i,k)
         LET xij = xx(i)-xx(j)
         LET yij = yy(i)-yy(j)
         LET rij = SQR(xij*xij+yij*yij)
         IF rij<rCutoff THEN
            LET f = force(rij,kind(i),kind(j))
            LET fxij = f*xij/rij
            LET fyij = f*yij/rij
            LET ffx(i) = ffx(i)+fxij
            LET ffy(i) = ffy(i)+fyij
            LET ffx(j) = ffx(j)-fxij
            LET ffy(j) = ffy(j)-fyij
         END IF
      NEXT k
   NEXT i
   FOR i=1 TO nMolec  ! boundary force
      LET ffx(i) = ffx(i)+boundaryForce(xx(i)+s)+boundaryForce(xx(i)-xMax-s)
      LET ffy(i) = ffy(i)+boundaryForce(yy(i)+s)+boundaryForce(yy(i)-yMax-s)
   NEXT i
END SUB

EXTERNAL FUNCTION force(r,ki,kj) !force(r) <-- forceTable - linear interporation
   LET ir = INT(r/hh)
   LET a = r - ir*hh
   LET force = ((hh-a)*forceTable(ki,kj,ir) + a*forceTable(ki,kj,ir+1))/hh
END FUNCTION

EXTERNAL FUNCTION boundaryForce(r)
   LET r6 = (3.418e-10/r)^6
   LET boundaryForce = (24.0*(0.5*1.711e-21)*r6*(2.0*r6-1.0)/r)
END FUNCTION

EXTERNAL SUB registerNearMolec
   LET rCut = rCutoff+20*2000*dt
   LET rcut2 = rCut*rCut
   FOR i=1 TO nMolec-1
      LET k = 1
      FOR j=i+1 TO nMolec
         LET r2 = (xx(i)-xx(j))*(xx(i)-xx(j))+(yy(i)-yy(j))*(yy(i)-yy(j))
         IF (r2<rcut2) THEN
            LET reg(i,k) = j
            LET k = k + 1
         END IF
      NEXT j
      LET reg(i,0) = k
   NEXT i
END SUB

EXTERNAL FUNCTION maxNearMolec
   LET mx = 0
   FOR i=1 TO nMolec-1
      IF mx<reg(i,0) THEN LET mx = reg(i,0)
   NEXT i
   LET maxNearMolec = mx-1
END FUNCTION

! ---------- utility

EXTERNAL FUNCTION systemTemprature
   LET kB = 1.38e-23 ! Boltzman's constant (J/K)
   LET ek= 0.0       !kinetic energy (J)
   FOR i=1 TO nMolec
      LET ek = ek + 0.5*mas(i)*(vx(i)^2+vy(i)^2)
   NEXT i
   LET systemTemprature = ek/(nMolec*kB) !2D: E/N=kT, 3D: E/N=(3/2)kT
END FUNCTION

EXTERNAL SUB ajustVelocity(temp)
   DECLARE EXTERNAL FUNCTION systemTemprature
   LET r = sqr(temp/systemTemprature)
   FOR i=1 TO nMolec
      LET vx(i) = r*vx(i)
      LET vy(i) = r*vy(i)
   NEXT i
END SUB

! ---------- draw particles

EXTERNAL SUB drawParticles(tempMode,contTemp,drawMode)
   DECLARE EXTERNAL FUNCTION systemTemprature,maxNearMolec
   DECLARE EXTERNAL sub realSpace,velocitySpace,plotBond
   SET DRAW MODE HIDDEN
   CLEAR
   IF drawMode=0 OR drawMode=1 THEN !--- 0:circle  1:circle+bond
      call plotBond(drawMode)
   ELSEIF drawMode=2 THEN
      call velocitySpace
   END IF
   !--- draw caption
   SET TEXT HEIGHT 10
   SET TEXT COLOR 1 ! black
   LET tmp$ = "adiabatic   constantTemp  "
   PLOT TEXT, AT  50, 70 ,USING "time =#####.## (ps)   temp =####.## (K)":sysTime*1E12,systemTemprature
   PLOT TEXT, AT  50, 55 ,USING molecSTR$(molecKind1)&","&molecSTR$(molecKind2)&"  N =####":nMolec
   PLOT TEXT, AT  50, 40 ,USING "tempMode = ## ":tempMode
   PLOT TEXT, AT 200, 40 :tmp$(tempMode*12+1:tempMode*12+12)
   PLOT TEXT, AT  50, 25 ,USING "controled Temperature =####.# (K)":contTemp
   PLOT TEXT, AT  50, 10 :"2-dimensional molecular dynamics"
   PLOT TEXT, AT  50,480 :"[esc] exit    [0] temp > value "
   PLOT TEXT, AT  50,460 :"[D] dispMode  [9] menu > select"
   SET DRAW MODE EXPLICIT
END SUB

EXTERNAL sub plotBond(drawMode)
   LET boxSize = 300
   LET mag = boxSize/xMax
   LET xp = 100
   LET yp = 100
   SET LINE COLOR 1 ! black : !--- box
   PLOT LINES: xp,yp; boxSize+xp,yp; boxSize+xp,boxSize+yp; xp,boxSize+yp; xp,yp
   SET TEXT HEIGHT 6
   SET TEXT COLOR 1 ! black
   PLOT TEXT, AT xp,boxSize+2+yp ,USING  "box size =##.# x ##.# (nm)":xMax*1e9,yMax*1e9
   FOR i=1 TO nMolec
      IF kind(i)=molecKind1 THEN SET LINE COLOR 11 ELSE SET LINE COLOR 12
      DRAW circle WITH SCALE(molecData(kind(i),3)/2.0*mag)*SHIFT(xx(i)*mag+xp,yy(i)*mag+yp)
   NEXT i
   IF drawMode=1 THEN
      FOR i=1 TO nMolec-1
         FOR k=1 TO reg(i,0)-1
            LET j = reg(i,k)
            LET r = SQR((xx(i)-xx(j))*(xx(i)-xx(j))+(yy(i)-yy(j))*(yy(i)-yy(j)))
            LET r0 = 0.5*(molecData(kind(i),3)+molecData(kind(j),3))
            IF r<1.2*r0 THEN
               LET yij = yy(i)-yy(j)
               IF yij=0.0 THEN LET yij = 1e-20
               LET th = 3.0*(ATN((xx(i)-xx(j))/yij)+PI/2.0)/PI
               LET col = 240 + INT((th-INT(th))*6)
               SET LINE COLOR col
               PLOT LINES: xx(i)*mag+xp,yy(i)*mag+yp;xx(j)*mag+xp,yy(j)*mag+yp
            END IF
         NEXT k
      NEXT i
   END IF
END sub

EXTERNAL sub velocitySpace
   LET boxSize = 300
   LET xp = 100
   LET yp = 100
   SET LINE COLOR 1 !black : axis
   PLOT LINES: xp,boxSize/2+yp; boxSize+xp,boxSize/2+yp !vx-axis
   PLOT LINES: boxSize/2+xp,yp; boxSize/2+xp,boxSize+yp !vy-axis
   SET TEXT HEIGHT 6
   SET TEXT COLOR 1 ! black
   PLOT TEXT, AT boxSize+xp,boxSize/2+yp: "vx"
   PLOT TEXT, AT boxSize+xp,boxSize/2-12+yp: "2000m/s"
   PLOT TEXT, AT boxSize/2-12+xp,boxSize+yp: "vy  2000m/s"
   PLOT TEXT, AT boxSize/2-8+xp,boxSize/2-10+yp: "0"
   PLOT TEXT, AT xp,boxSize+8+yp: "velocity space (vx,vy)"
   LET mag = boxSize/4000
   FOR i=1 TO nMolec
      IF kind(i)=molecKind1 THEN SET LINE COLOR 11 ELSE SET LINE COLOR 12
      DRAW circle WITH SCALE(5)*SHIFT(vx(i)*mag+boxSize/2+xp,vy(i)*mag+boxSize/2+yp)
   NEXT i
END sub
END MODULE
 

C/C++ 数学関数

 投稿者:しばっち  投稿日:2017年 6月18日(日)12時42分46秒
返信・引用
  C/C++で定義されている数学関数を定義してみました。

INPUT Z
LET L=LGAMMA(Z)
PRINT L;EXP(L)
END

EXTERNAL FUNCTION GAMMA(Z)
LET X$=PACKDBL$(Z)
IF Z>0 THEN
   LET GAMMA=GAMMA_(X$)
ELSE
   IF FP(Z)=0 THEN
      PRINT "定義域エラー"
      STOP
   END IF
   LET GAMMA=GAMMA_(X$)
END IF
FUNCTION GAMMA_(X$)
   ASSIGN ".\DLL\gamma_.dll","tgamma_",FPU
END FUNCTION
END FUNCTION

EXTERNAL FUNCTION LGAMMA(Z)
LET X$=PACKDBL$(Z)
LET LGAMMA=LGAMMA_(X$)

FUNCTION LGAMMA_(X$)
   ASSIGN ".\DLL\gamma_.dll","lgamma_",FPU
END FUNCTION
END FUNCTION

EXTERNAL FUNCTION DIGAMMA(Z)
LET X$=PACKDBL$(Z)
LET DIGAMMA=DIGAMMA_(X$)

FUNCTION DIGAMMA_(X$)
   ASSIGN ".\DLL\gamma_.dll","digamma_",FPU
END FUNCTION
END FUNCTION

EXTERNAL FUNCTION TRIGAMMA(Z)
LET X$=PACKDBL$(Z)
LET TRIGAMMA=TRIGAMMA_(X$)

FUNCTION TRIGAMMA_(X$)
   ASSIGN ".\DLL\gamma_.dll","trigamma_",FPU
END FUNCTION
END FUNCTION

EXTERNAL FUNCTION POLYGAMMA(N,Z)
LET X$=PACKDBL$(Z)
LET POLYGAMMA=POLYGAMMA_(INT(N),X$)

FUNCTION POLYGAMMA_(N,X$)
   ASSIGN ".\DLL\gamma_.dll","polygamma_",FPU
END FUNCTION
END FUNCTION
-------------------------------------------------------------------------------------------
                                        gamma_.cpp

#include <boost/math/special_functions/gamma.hpp>
#include <boost/math/special_functions/digamma.hpp>
#include <boost/math/special_functions/trigamma.hpp>
#include <boost/math/special_functions/polygamma.hpp>
using namespace boost::math;
using namespace std;

extern "C"  __declspec(dllexport)  double  tgamma_(double *a)
{
    return tgamma(*a);
}

extern "C"  __declspec(dllexport)  double  lgamma_(double *a)
{
    return lgamma(*a);
}

extern "C"  __declspec(dllexport)  double  digamma_(double *a)
{
    return digamma(*a);
}

extern "C"  __declspec(dllexport)  double  trigamma_(double *a)
{
    return trigamma(*a);
}

extern "C"  __declspec(dllexport)  double  polygamma_(int n,double *a)
{
    return polygamma(n,*a);
}



※ご注意

DLL側ではエラー処理をしていません。私自身が定義域をよく理解していませんので。
BASIC側で処理するか、範囲を超えないよう注意してください。



また、1000桁モードで使用できる関数も定義してみました。

OPTION ARITHMETIC DECIMAL_HIGH
LET X=1/5
CALL LATAN(1000,STR$(X),RESULT$)
LET S=16*VAL(RESULT$)
LET X=1/239
CALL LATAN(1000,STR$(X),RESULT$)
LET S=S-4*VAL(RESULT$)
PRINT S
PRINT PI
END

EXTERNAL  SUB LATAN(KETA,X$,RESULT$)
OPTION ARITHMETIC DECIMAL_HIGH
LET B$=REPEAT$(CHR$(0),KETA+100)
CALL ATAN1000(KETA,X$,B$)
FOR I=LEN(B$) TO 1 STEP -1
   IF B$(I:I)<="9" AND B$(I:I)>="0" THEN
      EXIT FOR
   END IF
NEXT I
LET RESULT$=B$(1:I)

SUB ATAN1000(KETA,X$,Y$)
   ASSIGN ".\DLL\atan1000.dll","atan1000"
END SUB
END SUB

EXTERNAL  SUB LACOTAN(KETA,X$,RESULT$)
OPTION ARITHMETIC DECIMAL_HIGH
LET B$=REPEAT$(CHR$(0),KETA+100)
CALL ACOTAN1000(KETA,X$,B$)
FOR I=LEN(B$) TO 1 STEP -1
   IF B$(I:I)<="9" AND B$(I:I)>="0" THEN
      EXIT FOR
   END IF
NEXT I
LET RESULT$=B$(1:I)

SUB ACOTAN1000(KETA,X$,Y$)
   ASSIGN ".\DLL\atan1000.dll","acotan1000"
END SUB
END SUB
-------------------------------------------------------------------------------------------
                                         atan1000.cpp

#include <cmath>
#include <string>
#include <boost/multiprecision/mpfr.hpp>

using namespace std;
using namespace boost::multiprecision;

extern "C"  __declspec(dllexport)  void atan1000(int keta,char *a,char *b)
{
    mpfr_float::default_precision(keta);
    mpfr_float x,y;
    string s;
    x.assign(a);
    y=atan(x);
    s=y.str(keta,true);
    strcpy(b,s.c_str());
}

extern "C"  __declspec(dllexport)  void acotan1000(int keta,char *a,char *b)
{
    mpfr_float::default_precision(keta);
    mpfr_float x,y;
    string s;
    x.assign(a);
    y=atan(1/x);
    s=y.str(keta,true);
    strcpy(b,s.c_str());
}


GMP互換ライブラリーを使用しています。

https://gmplib.org/

「mpir.dll」「mpfr.dll」をBASIC.EXEと同じフォルダーに入れてください。

http://mpir.org/
http://www.mpfr.org/


下記よりダウンロードしてください。(math.zip)

ダウンロードパス:shibacchi

http://fast-uploader.com/file/7053312569584/


なお、有効期限は本日より1ヶ月間となります。

P.S
VC++2017でコンパイルした Cファイルは問題なかったのですが、
C++ファイルはBASIC側で「DLLファイルがロードできない」等のエラーが出ました。
仕方なく諦めようかと思いましたが、VC++2010でコンパイルしてみたところ
問題なく動作しました。よってVC++2017とVC++2010でコンパイルしたDLLファイル
が混在しています。
 

Paract BASIC 0.9.1.1

 投稿者:白石和夫  投稿日:2017年 6月18日(日)08時53分30秒
返信・引用 編集済
  Paract BASIC 0.9.1.1で無用な再描画を削減しました。
描画に時間を要する比較的新しいバージョンのOS-Xでの実行で有効で,OS 10.9での所要時間がOS 10.6での所要時間との差が縮まっています。
なお、0.9.1Rev1でも同様ですが、OS-Xでは描画ループが回る時間が長くなるので、描画の実行途中が見れないことがあります。実行中での再描画が必要な場合は、SET DRAWMODE EXPLICIT を実行してください。

Intel Core 2 Duo (2.53GHz)
Mac OS 10.9.5

  Julia.BAS
  ParactBASIC 0.9.1.1で  3.09s

  Dragon.BAS
  ParactBASIC 0.9.1.1で 8.27s

結晶の分子動力学
  http://6317.teacup.com/basic/bbs/4308
ParactBASIC 0.9.1.1で 53.36s
 

Re: Paract BASIC 0.9.1

 投稿者:白石和夫  投稿日:2017年 6月17日(土)10時42分7秒
返信・引用
  MAC OS 10.9 での実行結果です。OS 10.6 より遅いです。

  Intel Core 2 Duo (2.53GHz)
  Mac OS 10.9.5

  Julia.BAS
  BASICAcc 1.0.4           で  5.40s
  ParactBASIC 0.9.1Rev1で  3.14s

  Dragon.BAS
  BASICAcc 1.0.4           で   3.41s
  ParactBASIC 0.9.1Rev1で 10.88s

結晶の分子動力学
  http://6317.teacup.com/basic/bbs/4308
BASICAcc 1.0.4           で122.60s
ParactBASIC 0.9.1Rev1で 84.36s
 

Paract BASIC 0.9.1

 投稿者:白石 和夫  投稿日:2017年 6月17日(土)10時19分25秒
返信・引用 編集済
  > No.4363[元記事へ]

Paract BASICをMAC OS-XとLinuxに対応させました。
描画と計算を別スレッドで実行するので 並行活動を使わないプログラムでもBASICAccより速くなることがあります。
ただし、 DRAW WITHで変形行列を計算させて自己相似図形を描くプログラムはBASICAccのほうが高速です。

Intel Core 2 Duo (2.53GHz)
Mac OS 10.6.8

Julia.BAS
BASICAcc 1.0.4           で  4.18s
ParactBASIC 0.9.1Rev1で  3.09s

Dragon.BAS
BASICAcc 1.0.4           で  4.34s
ParactBASIC 0.9.1Rev1で  6.22s

結晶の分子動力学
http://6317.teacup.com/basic/bbs/4308
BASICAcc 1.0.4           で 94.97s
ParactBASIC 0.9.1Rev1で 50.22s
 

cos(x)+ i sin(x)=ρ_n

 投稿者:たろさ  投稿日:2017年 6月13日(火)11時07分37秒
返信・引用 編集済
  (cos(x)+ i sin(x)=ρ_1)https://www.wolframalpha.com/

wolfram alpha変換(X) DATA の収束過程を観察するプログラム

訂正
※ 計算結果を見ると、実部の符号が全てマイナスになっています。

FOR n=0 TO 100
初期値が1になっていました。お詫びして訂正します。


(-.49999999999999  14.1347251417347)
              0.5 +14.13472514173469379045725(ρ_1)

(-.49999999999997  21.0220396387715)
              0.5 +21.02203963877155499262847(ρ_2)

(-.500000000000007 -25.0108575801457)
(-.499999999999899  25.0108575801456)
               0.5 +25.0108575801456887632137(ρ_3)

(-.500000000000205  30.4248761258594)
(-.499999999999996 -30.4248761258595)
(-.49999999999996  -30.4248761258596)
              0.5 + 30.4248761258595132103118(ρ_4)
--------------------------------------------------------

!e^(i*x)=ρ_n
!Σ(i*x)^n/FACT(n)=ρ_n

LET P=1 !step 1  (π周期)
LET x=2*(PI*P+COMPLEX(0.76771859766171526356,-1.32462990506382773366))!(cos(x)+ i sin(x)=ρ_1)https://www.wolframalpha.com/
!LET x=2*(PI*P+COMPLEX(-0.7677185976617152636,- 1.3246299050638277337))!(cos(x)- i sin(x)=ρ_1)
!LET x=2*(PI*P-COMPLEX( .767718597661715,1.32462990506383))!CONJ(+)ρ_1
!LET x=2*(PI*P+COMPLEX(0.7735081242682836818,- 1.5229270833490052893))!(cos(x)+ i sin(x)=ρ_2)
!LET x=2*(PI*P+COMPLEX(-0.7735081242682836818,- 1.5229270833490052893))!(cos(x)- i sin(x)=ρ_2)
!LET x=2*(PI*P-COMPLEX( .773508124268284,  1.52292708334901)) !CONJ(+)ρ_2
!LET x=2*(PI*P+COMPLEX(0.775403835822546139279,-1.609754910131780344447))!(cos(x)+ i sin(x)=ρ_3)
!LET x=2*(PI*P+COMPLEX(-0.7754038358225461393, - 1.6097549101317803444))!(cos(x)- i sin(x)=ρ_3)
!LET x=2*(PI*P-COMPLEX( .775403835822546,  1.60975491013178)  )!CONJ(+)ρ_3
!LET x=2*(PI*P+COMPLEX(0.7771819426816583383, - 1.7076977930242792999))!(cos(x)+ i sin(x)=ρ_4)
!LET x=2*(PI*P+COMPLEX(-0.7771819426816583383, -1.7076977930242792999))!(cos(x)- i sin(x)=ρ_4)
!LET x=2*(PI*P-COMPLEX( .777181942681658, 1.70769779302428)) !CONJ(+)ρ_4

!PRINT CONJ(COMPLEX(0.76771859766171526356,-1.32462990506382773366))
!PRINT CONJ(COMPLEX(0.7735081242682836818,- 1.5229270833490052893))
!PRINT CONJ(COMPLEX(0.775403835822546139279,-1.609754910131780344447))
!PRINT CONJ(COMPLEX(0.7771819426816583383, - 1.7076977930242792999))
LET i=COMPLEX(0,1)
FOR n=0 TO 100
   LET s=s+(i*x)^n/FACT(n)
   PRINT s
NEXT n
!PRINT s
PRINT
PRINT EXP(i*x)


END

---------------------------------------------------

!非自明な零点の変換
!e^(i*x)=ρ_n
!Σ(i*x)^n/FACT(n)=ρ_n
DATA 14.13472514173469379045725 !(ρ_1)
DATA 21.02203963877155499262847
DATA 25.0108575801456887632137
DATA 30.4248761258595132103118
DATA 32.9350615877391896906623
DATA 37.5861781588256712572177
DATA 40.9187190121474951873981
DATA 43.3270732809149995194961
DATA 48.0051508811671597279424
DATA 49.7738324776723021819167 !(ρ_10)

DIM A(10)
MAT READ A

LET P=1 !step 1  (π周期)
LET x=2*PI*P - COMPLEX(0,1)* LOG(COMPLEX(.5,A(1)))!(ρ_n)

LET i=COMPLEX(0,1)
FOR n=0 TO 80
   LET s=s+((i*x)^n/FACT(n))
   PRINT s
NEXT n
PRINT
PRINT EXP(i*x)
END

http://blogs.yahoo.co.jp/donald_stinger

 

Re: 計算結果が合いません。

 投稿者:たろさ  投稿日:2017年 6月12日(月)19時31分41秒
返信・引用
  > No.4378[元記事へ]

白石和夫さんへのお返事です。

> ご報告ありがとうございました。
> FACT関数は非負整数以外の入力に対しEXTYPE4000の例外を生成するのが正しいのですが,
> 2進モード,複素数モードでチェックが抜けていました。
> いずれ修正しますが,とりあえず引数が非負整数の範囲で使ってください。
> なお,これは,FACT(n)=PERM(n,n)で定義しているので,PERM関数のバグです。
> PERM(n,r)は,rが非負整数でないときEXTYPE4000になるべきでした。
>

お忙しいところ恐縮です。別件です。

----------------------------------------
!sin(x)+cos(x)=ρ_n ,sin(x)+cos(x)=zeta(n),
OPTION ARITHMETIC COMPLEX
SET BITMAP SIZE 640*2+1,640*2+1 ! x軸方向にデフォルトのn倍拡張
LET Wi=50
SET WINDOW -.5,2,-5,wi
SET COLOR MODE "REGULAR"
SET COLOR MIX(0) 0,0,0 !背景色黒
SET COLOR MIX(1) 0,0,1 !青
SET COLOR MIX(2) 1,0,0 !赤
SET COLOR MIX(3) 1,0,1 !赤紫
SET COLOR MIX(4) 0,1,0 !緑
SET COLOR MIX(5) 0,1,1 !水色
SET COLOR MIX(6) 1,1,0 !黄色
SET COLOR MIX(7) 1,1,1 !白
CLEAR !CLEAR文は,モード切替え前にSET COLOR MIX(0)で割り当てられた色で画面を塗りつぶす。
DRAW GRID(.5,10)
SET LINE COLOR 5
LET ZERO1=14.1435658457259942670362036372 !abs(ρ_1)
DRAW circle  WITH SCALE(ZERO1)*SHIFT(0,0)

LET ZERO2=21.0279849385071248276781391990
DRAW circle  WITH SCALE(ZERO2)*SHIFT(0,0)

LET ZERO3=25.0158549103229941605743395159
DRAW circle  WITH SCALE(ZERO3)*SHIFT(0,0)

LET ZERO4=30.4289843286609910063177407380
DRAW circle  WITH SCALE(ZERO4)*SHIFT(0,0)

LET ZERO5=32.9388567164704955328815401721
DRAW circle  WITH SCALE(ZERO5)*SHIFT(0,0)

LET A9=19

DATA 0.3751071671920218266   !ρ_n
DATA 0.3808358315180773568
DATA 0.3827206823919593268
DATA 0.3844917185921564522
DATA 0.3851159549678919465
DATA 0.3860527943271225999
DATA 0.3865933572251905934
DATA 0.3869323439242140375
DATA 0.3874937533934756098
DATA 0.3876785568691843285
DATA 0.392699081698724       !zeta(2)
DATA 0.11527703852715692264
DATA 0.0430686148913785
DATA 0.018822570327114491965
DATA 0.00874851187966394
DATA 0.0041922627745319701263
DATA 0.00204285704164369
DATA 0.0010052075320191288978
DATA 0.000497535187258165

DIM A(A9)
MAT READ A

DATA -1.4991587101393674264 !ρ_n
DATA -1.6967776724345700568
DATA -1.7834404027948773685
DATA -1.8812542252228525807
DATA -1.9208434596511019464
DATA -1.9868259982857292363
DATA -2.0292671172330924267
DATA -2.0578419991380527779
DATA -2.1090764852814033995
DATA -2.1271575682140434683
DATA -0.281861187900338      !zeta(2)
DATA 0
DATA 0
DATA 0
DATA 0
DATA 0
DATA 0
DATA 0
DATA 0

DIM B(A9)
MAT READ B
SET POINT STYLE 7
LET i=COMPLEX(0,1)
FOR n=1 TO A9
   FOR nn=1 TO 4
      IF n<11 THEN
         LET x=2*(PI*nn - COMPLEX(A(n) ,B(n))) !ρ_n
      ELSE
         LET x=2*(PI*nn + COMPLEX(A(n) ,B(n))) !zeta(n)
      END IF
      LET z1=(1/2 + i/2)*EXP(-i*x) + (1/2 - i/2)*EXP(i*x)
      LET s1=(1/2 + i/2)*EXP(-i*x)
      LET s2=(1/2 - i/2)*EXP(i*x)
      SET POINT COLOR 2
      PLOT POINTS: RE(s1)+RE(s2) ,IM(s1)+IM(s2)
      PRINT ROUND(n,2);":";z1;RE(s1)+RE(s2);IM(s1)+IM(s2)
   NEXT nn
NEXT n
END
------------------------------------------


LET x=2*(PI*nn - COMPLEX(A(n) ,B(n))) !ρ_n

LET x=2*(PI*nn + COMPLEX(A(n) ,B(n))) !zeta(n)


X  DATA

sin(x)+cos(x)=zeta(n)
sin(x)+cos(x)=ρ_n
https://www.wolframalpha.com/


LET x=2*(PI*nn (-) COMPLEX(A(n) ,B(n))) !ρ_n

LET x=2*(PI*nn (+) COMPLEX(A(n) ,B(n))) !zeta(n)

この ± のギャップ 高精度計算サイトでも経験しました。



http://blogs.yahoo.co.jp/donald_stinger

 

Re: 計算結果が合いません。

 投稿者:白石和夫  投稿日:2017年 6月12日(月)18時18分15秒
返信・引用 編集済
  > No.4377[元記事へ]

ご報告ありがとうございました。
FACT関数は非負整数以外の入力に対しEXTYPE4000の例外を生成するのが正しいのですが,
2進モード,複素数モードでチェックが抜けていました。
いずれ修正しますが,とりあえず引数が非負整数の範囲で使ってください。
なお,これは,FACT(n)=PERM(n,n)で定義しているので,PERM関数のバグです。
PERM(n,r)は,rが非負整数でないときEXTYPE4000になるべきでした。



 

計算結果が合いません。

 投稿者:たろさ  投稿日:2017年 6月12日(月)12時30分41秒
返信・引用
  OPTION ARITHMETIC COMPLEX

FOR n=2 TO 10 STEP 0.1
   LET x=LOG(FACT(n-1))
   PRINT ROUND(n,2);x
NEXT n
END
---------------------------------
BASIC Acc と 十進BASIC

高精度計算サイトhttp://keisan.casio.jp/calculator
https://www.wolframalpha.com/
ln(gamma(x));
ln((x-1)!)


http://blogs.yahoo.co.jp/donald_stinger

 

Re: トーラス

 投稿者:白石和夫  投稿日:2017年 6月 4日(日)18時17分18秒
返信・引用
  > No.4375[元記事へ]

白石和夫さんへのお返事です。

> >
> > 注釈外すと正常に見えますが、バー(2)を動かすとその値がバー(12)にも表示されます。
> >
> > !LOCATE VALUE NOWAIT(2),RANGE -5 TO 5,AT 0:B
> > LOCATE VALUE NOWAIT(11),RANGE -5 TO 5,AT 0:K
> > LOCATE VALUE NOWAIT(12),RANGE -5 TO 5,AT 0:L
> > LOCATE VALUE NOWAIT(13),RANGE -5 TO 5,AT 0:M
> > DO
> > !LOCATE VALUE NOWAIT(2): B
> >    LOCATE VALUE NOWAIT(11): K
> >    LOCATE VALUE NOWAIT(12): L
> >    LOCATE VALUE NOWAIT(13): M
> > LOOP
> > END
> >
>
> 現象を確認しました。
> どこか間違えている気がします。
>
間違えている箇所を特定しました。
順次,修正版を用意します。

 

Re: トーラス

 投稿者:白石和夫  投稿日:2017年 6月 4日(日)18時04分50秒
返信・引用
  > No.4374[元記事へ]

>
> 注釈外すと正常に見えますが、バー(2)を動かすとその値がバー(12)にも表示されます。
>
> !LOCATE VALUE NOWAIT(2),RANGE -5 TO 5,AT 0:B
> LOCATE VALUE NOWAIT(11),RANGE -5 TO 5,AT 0:K
> LOCATE VALUE NOWAIT(12),RANGE -5 TO 5,AT 0:L
> LOCATE VALUE NOWAIT(13),RANGE -5 TO 5,AT 0:M
> DO
> !LOCATE VALUE NOWAIT(2): B
>    LOCATE VALUE NOWAIT(11): K
>    LOCATE VALUE NOWAIT(12): L
>    LOCATE VALUE NOWAIT(13): M
> LOOP
> END
>

現象を確認しました。
どこか間違えている気がします。
 

Re: トーラス

 投稿者:しばっち  投稿日:2017年 6月 4日(日)14時37分1秒
返信・引用
  > No.4373[元記事へ]

たろささんへのお返事です。

> 毎度、力不足で申し訳ありませんが、トーラス型の半球の継ぎ手が出てます。

浮動小数による周期のずれか、誤差が原因ではないかと思います。それ以上は分かりません。



> 出来たら、二重螺旋も出来たらと、思います。

方程式を訂正すれば可能だとは思いますが...

http://izumi-math.jp/sanae/MathCurves/MathCurves.htm



> 家のパソコンでは、LOCATE VALUE NOWAIT(12)が調子悪い。別プログラムで確かめましたが
> 特に問題ナシです。
> Lazarus 1.6.4 (FPC 3.0.2) BASIC Accelerator Ver. 1.0.3 を使用しました。

Windows Lazarus版 Basic Ver6.6.3.0 でもLOCATE VALUE NOWAIT(12)のおかしなところ確認できました。

注釈外すと正常に見えますが、バー(2)を動かすとその値がバー(12)にも表示されます。

!LOCATE VALUE NOWAIT(2),RANGE -5 TO 5,AT 0:B
LOCATE VALUE NOWAIT(11),RANGE -5 TO 5,AT 0:K
LOCATE VALUE NOWAIT(12),RANGE -5 TO 5,AT 0:L
LOCATE VALUE NOWAIT(13),RANGE -5 TO 5,AT 0:M
DO
!LOCATE VALUE NOWAIT(2): B
   LOCATE VALUE NOWAIT(11): K
   LOCATE VALUE NOWAIT(12): L
   LOCATE VALUE NOWAIT(13): M
LOOP
END


 

Re: トーラス

 投稿者:たろさ  投稿日:2017年 6月 4日(日)09時25分47秒
返信・引用 編集済
  > No.4372[元記事へ]

たろささんへのお返事です。

連続で失礼します。

トーラス型から球体や螺旋型へ変形出来ると思い挑戦しました。

!LOCATE VALUE スライドバーを12本使用しています
!lazaus版BASIC.exe ver 6.6.3.0 又はParact Basic,Basic Accを使用してください
OPTION ARITHMETIC COMPLEX
!OPTION ARITHMETIC NATIVE
OPTION ANGLE DEGREES
RANDOMIZE
!ASK MAX VALUE DEVICE nv
!PRINT nv
LET zi=COMPLEX(0,1)
LET ZTH=0          ! z軸のまわりの回転角
LET XTH=0          ! x軸のまわりの回転角初期値
LET YTH=0          ! y軸のまわりの回転角初期値
LET LMIN=1E+10
LET LMAX=-1E+10
LET NN=40 !'分割数
LET MM=40
DIM M(4,4),POINT(4),ROTX(4,4),ROTY(4,4)
DIM XX(0 TO NN,0 TO MM),YY(0 TO NN,0 TO MM),ZZ(0 TO NN,0 TO MM)
MAT M=ROTATE(ZTH)
SET POINT STYLE 1
LET XDT=RND-.5
LET YDT=RND-.5
LOCATE VALUE NOWAIT(1),RANGE 0 TO 2,AT 1 : SCALE
LOCATE VALUE NOWAIT(2),RANGE -3 TO 3,AT 1 : SPEED
LOCATE VALUE NOWAIT(3),RANGE 0 TO 20,AT 0.5 : K
LOCATE VALUE NOWAIT(4),RANGE 0 TO 20,AT 1.5 : RR
LOCATE VALUE NOWAIT(5),RANGE 0 TO 20,AT 1.5 : R0
LOCATE VALUE NOWAIT(6),RANGE 0 TO 20,AT 5.4 : R1
LOCATE VALUE NOWAIT(7),RANGE 0 TO 20,AT 1.5 : R2
LOCATE VALUE NOWAIT(8),RANGE -100 TO 100,AT 0:XMOVE
LOCATE VALUE NOWAIT(9),RANGE -100 TO 100,AT 0:YMOVE
LOCATE VALUE NOWAIT(10),RANGE -100 TO 100,AT 0:ZMOVE
LOCATE VALUE NOWAIT(11),RANGE 1 TO 20,AT 3.58 : Wa
LOCATE VALUE NOWAIT(13),RANGE 1 TO 20,AT 3.58 : Wb
DO
   LOCATE VALUE NOWAIT(1): SCALE
   LOCATE VALUE NOWAIT(2): SPEED
   LOCATE VALUE NOWAIT(3): K
   LOCATE VALUE NOWAIT(4): RR
   LOCATE VALUE NOWAIT(5): R0
   LOCATE VALUE NOWAIT(6): R1
   LOCATE VALUE NOWAIT(7): R2
   LOCATE VALUE NOWAIT(8): XMOVE
   LOCATE VALUE NOWAIT(9): YMOVE
   LOCATE VALUE NOWAIT(10): ZMOVE
   LOCATE VALUE NOWAIT(11): Wa
   LOCATE VALUE NOWAIT(13): Wb
   LET LL=0
   FOR J=0 TO MM
      FOR I=0 TO NN
         LET LL=LL+0.00383
         LET ALPHA=I*360/NN     !トーラス
         LET BETA=J*360/MM
         LET  XX(I,J)=(R1+R0*wa*COS(ALPHA))*(R2+RR*COS(K*BETA))*RE(EXP(LL*zi))!*COS(BETA)
         LET  ZZ(I,J)=(R1+R0*wb*COS(ALPHA))*(R2+RR*COS(K*BETA))*IM(EXP(LL*zi))!*SIN(BETA)
         LET  YY(I,J)=R0*SIN(ALPHA)*(R2+RR*COS(K*BETA))
      NEXT I
   NEXT  J
   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
   SET DRAW MODE HIDDEN
   CLEAR
   FOR J=0 TO MM-1
      FOR I=0 TO NN-1
         CALL PLOT(XX(I,J),YY(I,J),ZZ(I,J))
         CALL PLOT(XX(I,J+1),YY(I,J+1),ZZ(I,J+1))
         CALL PLOT(XX(I+1,J+1),YY(I+1,J+1),ZZ(I+1,J+1))
         CALL PLOT(XX(I+1,J),YY(I+1,J),ZZ(I+1,J))
         CALL PLOT(XX(I,J),YY(I,J),ZZ(I,J))
         PLOT LINES
      NEXT   I
   NEXT  J
   IF FL=0 THEN
      SET WINDOW -LMAX*1.2,LMAX*1.2,-LMAX*1.2,LMAX*1.2
      LET WW=LMAX*2.4
   END IF
   LET FL=1
   SET DRAW MODE EXPLICIT
   MOUSE POLL X,Y,L,R
   IF R<>0 THEN EXIT DO
   LET XTH=0
   LET YTH=0
   IF L<>0 THEN
      DO WHILE L<>0
         MOUSE POLL X,Y,L,R
      LOOP
      LET XDT=-(Y-Y0)/WW*5
      LET YDT= (X-X0)/WW*5
      LET XDT=MAX(-5,MIN(5,XDT))
      LET YDT=MAX(-5,MIN(5,YDT))
   ELSE
      LET XTH=XTH+XDT*SPEED
      LET YTH=YTH+YDT*SPEED
   END IF
   LET X0=X
   LET Y0=Y
LOOP

SUB PLOT(X,Y,Z)

   LET POINT(1)=X+XMOVE
   LET POINT(2)=Y+YMOVE
   LET POINT(3)=Z+ZMOVE
   MAT POINT=POINT*M
   IF FL=0 THEN
      LET LMIN=MIN(LMIN,POINT(1))
      LET LMAX=MAX(LMAX,POINT(1))
      LET LMIN=MIN(LMIN,POINT(2))
      LET LMAX=MAX(LMAX,POINT(2))
   ELSE
      SET LINE COLOR 2
      PLOT LINES:POINT(1)*SCALE,POINT(2)*SCALE;
   END IF
END SUB
END

--------------------------------------------

初期値はオオム貝型


螺旋型
k 0.52
RR 1.54
R0 1.54
R1 5.54
R2 1.54

ゼータ関数の桃型
k 1.3
RR 1.5
R0 1.5
R1 5.4
R2 1.5

フェルマー螺旋型
k 0.24
RR 3.52
R0 3.52
R1 0
R2 0

ハート型
k 1.02
RR 1.54
R0 1.54
R1 5.54
R2 1.54

毎度、力不足で申し訳ありませんが、トーラス型の半球の継ぎ手が出てます。

出来たら、二重螺旋も出来たらと、思います。

ご協力をお願いします。


家のパソコンでは、LOCATE VALUE NOWAIT(12)が調子悪い。別プログラムで確かめましたが
特に問題ナシです。
Lazarus 1.6.4 (FPC 3.0.2) BASIC Accelerator Ver. 1.0.3 を使用しました。
LACATE VALUE スライドバー個数の上限20本を確認です。

http://blogs.yahoo.co.jp/donald_stinger

 

Re: トーラス

 投稿者:たろさ  投稿日:2017年 6月 3日(土)21時40分1秒
返信・引用
  > No.4371[元記事へ]

しばっちさんへのお返事です。

> たろささんへのお返事です。
>
> > トーラスを球体付近から、平面付近にに、変形出来たら、と思いました。
>
> LOCATE VALUE スライドバーを10本使用しています
> lazaus版BASIC.exe ver 6.6.3.0 又はParact Basic,Basic Accを使用してください
>

有難うございました。

Paract Basic Ver.902  動作確認

ゼータ関数の非自明な零点グラフと似た図形が取れました。

OPTION ARITHMETIC COMPLEX
LET wi=3
SET WINDOW -wi+1,wi+1,-wi,wi!-2.5,2.5,-2.5,2.5 !x,y
SET COLOR MODE "REGULAR"
SET COLOR MIX(0) 0,0,0 !背景色黒
SET COLOR MIX(1) 0,0,1 !青
SET COLOR MIX(2) 1,0,0 !赤
SET COLOR MIX(3) 0,1,0 !緑
SET COLOR MIX(4) 1,0,1 !赤紫
SET COLOR MIX(5) 0,1,1 !水色
SET COLOR MIX(6) 1,1,0 !黄色
SET COLOR MIX(7) 1,1,1 !白
CLEAR !CLEAR文は,モード切替え前にSET COLOR MIX(0)で割り当てられた色で画面を塗りつぶす。

DRAW grid
SET POINT STYLE 7
LET c=0
FOR x=30.42487 TO 40.91871 STEP 0.01 !ρ_4~ρ_7
   IF x=1 THEN GOTO 100
   LET m=Zeta(COMPLEX(1/2,X))
   LET R=RE(m)
   LET I=IM(m)
   ! PRINT m
   LET c=c+1
   LET Ca=MOD(c,3)+1
   SET POINT COLOR Ca
   PLOT POINTS:r,i
   PLOT
100 NEXT x
END

EXTERNAL  FUNCTION ZETA(S) !'ゼータ関数
    OPTION ARITHMETIC COMPLEX
    FOR M=1 TO 300
       LET SS=0
       FOR J=1 TO M
          LET SS=SS+(-1)^(J-1)*COMB(M-1,J-1)*J^(-S)
       NEXT J
       LET SUM=SUM+SS*2^(-M)
       IF ABS(SUM-S0)<1E-14 THEN
          LET ZETA=SUM/(1-2^(1-S))
          EXIT FUNCTION
       END IF
       LET S0=SUM
    NEXT M
    PRINT "収束エラー"
    STOP
END FUNCTION

http://blogs.yahoo.co.jp/donald_stinger

 

Re: トーラス

 投稿者:しばっち  投稿日:2017年 6月 3日(土)18時43分17秒
返信・引用
  > No.4370[元記事へ]

たろささんへのお返事です。

> トーラスを球体付近から、平面付近にに、変形出来たら、と思いました。

LOCATE VALUE スライドバーを10本使用しています
lazaus版BASIC.exe ver 6.6.3.0 又はParact Basic,Basic Accを使用してください

OPTION ARITHMETIC NATIVE
OPTION ANGLE DEGREES
RANDOMIZE
LET ZTH=0          ! z軸のまわりの回転角
LET XTH=0          ! x軸のまわりの回転角初期値
LET YTH=0          ! y軸のまわりの回転角初期値
LET LMIN=1E+10
LET LMAX=-1E+10
LET NN=40 !'分割数
LET MM=40
DIM M(4,4),POINT(4),ROTX(4,4),ROTY(4,4)
DIM XX(0 TO NN,0 TO MM),YY(0 TO NN,0 TO MM),ZZ(0 TO NN,0 TO MM)
MAT M=ROTATE(ZTH)
SET POINT STYLE 1
LET XDT=RND-.5
LET YDT=RND-.5
LOCATE VALUE NOWAIT(1),RANGE 0 TO 2,AT 1 : SCALE
LOCATE VALUE NOWAIT(2),RANGE -3 TO 3,AT 1 : SPEED
LOCATE VALUE NOWAIT(3),RANGE 0 TO 20,AT 1 : K
LOCATE VALUE NOWAIT(4),RANGE 0 TO 20,AT 0 : RR
LOCATE VALUE NOWAIT(5),RANGE 0 TO 20,AT 5 : R0
LOCATE VALUE NOWAIT(6),RANGE 0 TO 20,AT 10 : R1
LOCATE VALUE NOWAIT(7),RANGE 0 TO 20,AT 1 : R2
LOCATE VALUE NOWAIT(8),RANGE -100 TO 100,AT 0:XMOVE
LOCATE VALUE NOWAIT(9),RANGE -100 TO 100,AT 0:YMOVE
LOCATE VALUE NOWAIT(10),RANGE -100 TO 100,AT 0:ZMOVE
DO
   LOCATE VALUE NOWAIT(1): SCALE
   LOCATE VALUE NOWAIT(2): SPEED
   LOCATE VALUE NOWAIT(3): K
   LOCATE VALUE NOWAIT(4): RR
   LOCATE VALUE NOWAIT(5): R0
   LOCATE VALUE NOWAIT(6): R1
   LOCATE VALUE NOWAIT(7): R2
   LOCATE VALUE NOWAIT(8): XMOVE
   LOCATE VALUE NOWAIT(9): YMOVE
   LOCATE VALUE NOWAIT(10): ZMOVE
   FOR J=0 TO MM
      FOR I=0 TO NN
         LET ALPHA=I*360/NN
         LET BETA=J*360/MM
         LET  XX(I,J)=(R1+R0*COS(ALPHA))*(R2+RR*COS(K*BETA))*COS(BETA)
         LET  ZZ(I,J)=(R1+R0*COS(ALPHA))*(R2+RR*COS(K*BETA))*SIN(BETA)
         LET  YY(I,J)=R0*SIN(ALPHA)*(R2+RR*COS(K*BETA))
      NEXT I
   NEXT  J
   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
   SET DRAW MODE HIDDEN
   CLEAR
   FOR J=0 TO MM-1
      FOR I=0 TO NN-1
         CALL PLOT(XX(I,J),YY(I,J),ZZ(I,J))
         CALL PLOT(XX(I,J+1),YY(I,J+1),ZZ(I,J+1))
         CALL PLOT(XX(I+1,J+1),YY(I+1,J+1),ZZ(I+1,J+1))
         CALL PLOT(XX(I+1,J),YY(I+1,J),ZZ(I+1,J))
         CALL PLOT(XX(I,J),YY(I,J),ZZ(I,J))
         PLOT LINES
      NEXT   I
   NEXT  J
   IF FL=0 THEN
      SET WINDOW -LMAX*1.2,LMAX*1.2,-LMAX*1.2,LMAX*1.2
      LET WW=LMAX*2.4
   END IF
   LET FL=1
   SET DRAW MODE EXPLICIT
   MOUSE POLL X,Y,L,R
   IF R<>0 THEN EXIT DO
   LET XTH=0
   LET YTH=0
   IF L<>0 THEN
      DO WHILE L<>0
         MOUSE POLL X,Y,L,R
      LOOP
      LET XDT=-(Y-Y0)/WW*5
      LET YDT= (X-X0)/WW*5
      LET XDT=MAX(-5,MIN(5,XDT))
      LET YDT=MAX(-5,MIN(5,YDT))
   ELSE
      LET XTH=XTH+XDT*SPEED
      LET YTH=YTH+YDT*SPEED
   END IF
   LET X0=X
   LET Y0=Y
LOOP

SUB PLOT(X,Y,Z)
   LET POINT(1)=X+XMOVE
   LET POINT(2)=Y+YMOVE
   LET POINT(3)=Z+ZMOVE
   MAT POINT=POINT*M
   IF FL=0 THEN
      LET LMIN=MIN(LMIN,POINT(1))
      LET LMAX=MAX(LMAX,POINT(1))
      LET LMIN=MIN(LMIN,POINT(2))
      LET LMAX=MAX(LMAX,POINT(2))
   ELSE
      PLOT LINES:POINT(1)*SCALE,POINT(2)*SCALE;
   END IF
END SUB
END
 

Re: トーラス

 投稿者:たろさ  投稿日:2017年 6月 2日(金)01時24分24秒
返信・引用
  > No.3135[元記事へ]

しばっちさんへのお返事です。

リーマン予想の探究

私は素人なので、普通の方法ではありません。

オイラーの公式(e^(i*pi)=-1)と(sin(pi)+cos(pi)=-1)から

x=1.080929170560724870066442798932725233917386163162;
s=2;
sqrt(sin(x)*s+cos(x)*s);
zeta(2);
sqrt((sin(x)*s)^2+(cos(x)*s)^2);
sin(x)*s;
cos(x)*s;

1.6449340668482264364724151666460251892189499012068=zeta(2)
1.6449340668482264364724151666460251892189499012068
2
1.7647907410794978581738952982401046226046426494767
0.94101734319834762061611394311281513433223473032

ゼータ関数を座標値に変換しました。
ゼータ関数の円周上の点 グラフv3
https://blogs.yahoo.co.jp/donald_stinger/15645572.html

------------------------------------------------------

(sin(x)*s+cos(x)*s)から、非自明な零点を求めて見たいと思いました。

n*pi/4    MOD(n,4)=3 の時 0

2/(pi*euler)=1.10291492604619838  ,(euler=0.577215664901532861)

等が原因で実験失敗


実験の副産物の報告です。

リーマン予想とは?  https://blogs.yahoo.co.jp/donald_stinger/15656946.html
-------------------------------------
実部は変化しない。
----------------------
y=sin(pi*x)+x*(1-x)i;
Re(y);
Im(y);
----------------------
-5,0.5,101
----------------------
y=sin(pi*x)-x*(1-x)i;
Re(y);
Im(y);
----------------------
-5,0.5,101
----------------------

実験プログラムです。
----------------------------------------
!v1
OPTION ARITHMETIC COMPLEX
SET BITMAP SIZE 640*1+1,640*1+1 ! x軸方向にデフォルトのn倍拡張
LET Wi=10
SET WINDOW -wi,wi,-wi,wi
DRAW GRID(1,1)
SET POINT STYLE 7
LET s=COMPLEX(0,1)
!LET s=COMPLEX(0.5,14.13472514173469379)
!LET s=2
LET z=0
FOR n=1 TO 100 STEP 0.1 !半径
   LET Z=Z+1/n*s
   LET z1=RE(z)
   LET z2=IM(z)
   SET POINT COLOR 2
   PLOT POINTS:z1+0.5,z2
NEXT n
PRINT z
END
-----------------------------------------------
!v2
OPTION ARITHMETIC COMPLEX
SET BITMAP SIZE 640*1+1,640*1+1 ! x軸方向にデフォルトのn倍拡張
LET Wi=10
SET WINDOW -wi,wi,-wi,wi
DRAW GRID(1,1)
SET POINT STYLE 7
!LET s=COMPLEX(0,1)
LET s=COMPLEX(.5,14.13472514173469379)
!LET s=2
LET z=0
FOR n=1 TO 100 STEP 0.1 !半径
   LET Z=Z+1/n^s
   LET z1=RE(z)
   LET z2=IM(z)
   SET POINT COLOR 2
   PLOT POINTS:z1,z2
NEXT n
PRINT z
END
------------------------------------------------
!v3
OPTION ARITHMETIC COMPLEX
SET BITMAP SIZE 640*1+1,640*1+1 ! x軸方向にデフォルトのn倍拡張
LET Wi=15
SET WINDOW -1,1,13,15
DRAW GRID(.1,.1)
SET POINT STYLE 7
!LET s=COMPLEX(0,1)
LET s=COMPLEX(.5,14.13472514173469379)
!LET s=2
LET z=0
FOR n=1 TO 100 STEP 0.1 !半径
   LET Z=Z+1/n*s
   LET z1=RE(z)
   LET z2=IM(z)
   SET POINT COLOR 2
   PLOT POINTS:z1,z2
NEXT n
PRINT z
END
------------------------------------------------
!v4
OPTION ARITHMETIC COMPLEX
SET BITMAP SIZE 640*1+1,640*1+1 ! x軸方向にデフォルトのn倍拡張
LET Wi=10
SET WINDOW -wi,wi,-wi,wi
DRAW GRID(1,1)
SET POINT STYLE 7

LET s=COMPLEX(0.5,14.13472514173469379)
!LET s=2
LET z=0
FOR n=1 TO 100 STEP 0.1 !半径
   LET Z=Z+1/n^s
   LET z1=RE(z)
   LET z2=IM(z)
   SET POINT COLOR 2
   PLOT POINTS:z2-z1,z1-z2
   SET POINT COLOR 4
   PLOT POINTS:z1-z2,z2-z1
NEXT n
PRINT z
END
-------------------------------------------

!v2 螺旋のグラフを見ると、カオスと似ている。

LET s=COMPLEX(z)  変数にすると螺旋の形状が変化する。

カオスの例に倣うと、物理的な現象では、と、予想すると

非自明な零点は、位相が逆転する点と予想される。

証明には物理学の知識が必要なのかも ?

次の実験の場をトーラスにして、実験中です。

現状かなり複雑です。



トーラス  http://6317.teacup.com/basic/bbs/3135

回転体   http://6317.teacup.com/basic/bbs/1621
Re: 回転体 http://6317.teacup.com/basic/bbs/1622

大変便利です。3Dプリンターの時代にフィットしてます。


Re: 平行投影 http://6317.teacup.com/basic/bbs/4365

トーラスを組み込んで、みましたが、スライドバーのためか?  失敗しています。

私には、3Dは、難しい。x,y,z のz軸が わかりません。

トーラスを球体付近から、平面付近にに、変形出来たら、と思いました。

http://blogs.yahoo.co.jp/donald_stinger

 

レンタル掲示板
/150