2.3 陽解法 heat1d-e.bas

もっとも簡単な陽解法のプログラムを以下に紹介する (この短さが BASIC の嬉しいところである)。

ちなみにグラフィックスがらみの命令は SET WINDOW, PLOT LINES, PLOT TEXT の 3 つのみ (最後の1つはなくてもグラフは描ける)。

REM heat1d-e.bas
REM 陽解法による１次元熱方程式(同次ディリクレ境界条件)
REM u_t(x,t)=u_{xx}(x,t) in (0,1)×(0,∞)
REM u(0,t)=u(1,t)=0, t∈(0,∞)
REM u(x,0)=f(x), x∈[0,1]
REM --------------------------------------------------
DECLARE EXTERNAL FUNCTION f
REM データ入力
INPUT PROMPT "分割数: ": N
INPUT PROMPT "λ: ": lambda
INPUT PROMPT "最終時刻: ": Tmax
DIM u(0 TO N),newu(0 TO N)
SET WINDOW -0.2,1.2,-0.2,1.2
PLOT TEXT ,AT 0.4,0.8: "熱方程式"
PLOT TEXT ,AT 0.4,0.7, USING "N=###,λ=##.####":N,lambda
LET h=1/N
LET tau=lambda*h*h
LET kmax=Tmax/tau
FOR i=0 TO N
   LET u(i)=f(i*h)
   PRINT u(i)
NEXT i
FOR i=0 TO N
   PLOT LINES: i*h,u(i);
NEXT i
FOR k=1 TO kmax
   FOR i=1 TO N-1
      LET newu(i)=(1-2*lambda)*u(i)+lambda*(u(i+1)+u(i-1))
   NEXT i
   FOR i=1 TO N-1
      LET u(i)=newu(i)
   NEXT i
   LET u(0)=0.0
   LET u(N)=0.0
   FOR i=0 TO N
      PLOT LINES: i*h,u(i);
   NEXT i
NEXT k
END
REM --------------------------------------------------
REM 初期データ
EXTERNAL FUNCTION f(x)
IF (x < 0.5) THEN
   LET  f=x
ELSE
   LET  f=1-x
END IF
END FUNCTION