2.4 θ法 heat1d-i.bas

REM heat1d-i.bas
REM 陰解法による１次元熱方程式 (同次ディリクレ境界条件)
REM u_t(x,t)=u_{xx}(x,t), (x,t)∈(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
DECLARE EXTERNAL SUB trilu,trisol
INPUT PROMPT "分割数: ": N
IF N <= 0 THEN
   LET  N=20
   LET  lambda=0.5
   LET  theta=0.5
   LET  Tmax=1
   LET  dt=0.01
ELSE
   INPUT PROMPT "λ: ": lambda
   INPUT PROMPT "θ: ": theta
   INPUT PROMPT "最終時刻: ": Tmax
   INPUT PROMPT "Δt: ": dt
END IF
DIM u(0 TO N),ff(N-1),al(N-1),ad(N-1),au(N-1)
FOR i=1 TO N-1
   LET  al(i)=-theta*lambda
   LET  ad(i)=1+2*theta*lambda
   LET  au(i)=-theta*lambda
NEXT i
CALL trilu(N-1,al,ad,au)
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,theta
LET  h=1/N
LET  tau=lambda*h*h
LET  SKIP=INT(dt/tau)
LET  nmax=Tmax/tau
FOR i=0 TO N
   LET u(i)=f(i*h)
NEXT i
FOR i=0 TO N
   PLOT LINES: i*h,u(i);
NEXT i
LET  c3=1-2*(1-theta)*lambda
LET  c4=(1-theta)*lambda
FOR k=1 TO nmax
   FOR i=1 TO N-1
      LET ff(i)=c3*u(i)+c4*(u(i+1)+u(i-1))
   NEXT i
   CALL trisol(N-1,al,ad,au,ff)
   FOR i=1 TO N-1
      LET  u(i)=ff(i)
   NEXT i
   LET u(0)=0.0
   LET u(N)=0.0
   IF MOD(k,SKIP)=0 THEN
      FOR i=0 TO N
         PLOT LINES: i*h,u(i);
      NEXT i
   END if
NEXT k
END
REM **************** 初期データ *****************
EXTERNAL FUNCTION f(x)
IF (x < 0.5) THEN
   LET  f=x
ELSE
   LET  f=1-x
END IF
END FUNCTION
REM *************** LU分解 **********************
EXTERNAL SUB trilu(n,al(),ad(),au())
FOR i=1 TO n-1
   LET al(i+1) = al(i+1) / ad(i)
   LET ad(i+1) = ad(i+1) - au(i) * al(i+1)
NEXT i
END SUB
REM *********************************************
EXTERNAL SUB trisol(n,al(),ad(),au(),b())
FOR i=1 TO n-1
   LET  b(i+1) = b(i+1) - b(i) * al(i+1)
NEXT i
LET  b(n) = b(n) / ad(n)
FOR i=n-1 TO 1 STEP -1
   LET  b(i) = (b(i) - au(i) * b(i+1)) / ad(i)
NEXT i
END SUB