吸収　ホームページプログラミング　

VBScript部はルンゲ・クッタ法による連立常微分方程式の解法プログラム。absorp.html

<SCRIPT language=VBScript><!--
Option Explicit
Dim vbCrLf : vbCrLf = chr(13) & chr(10)
Dim Y(), Y0(), YW(), YS(), YQ(), YD()
Dim NM,I,YA,yp2,y1,y2,CL,YB,KYA,KXA,G,L,x1,M,D,YH0,YE0,YH1,X0,NK,NZ,YH,X,YER,E1,EN,YE1,YDD,DYDX,y20,Z
Dim GZ(30),GY(30),GYI(30),GX(30),GXI(30)
Sub A_onClick
txtsolution.value=""
document.chartApplet.clear
chartApplet.setColor 255,255,200
chartApplet.fillRect 30,10,300,300
chartApplet.setColor 0,0,0
chartApplet.setFontSize 12
chartApplet.drawLine 30, 10, 330, 10
chartApplet.drawLine 30, 310, 330, 310
chartApplet.drawLine 330, 10, 330, 310
chartApplet.drawLine 30, 10, 30, 310
chartApplet.drawLine 80, 300, 80, 310
chartApplet.drawLine 130, 300, 130, 310
chartApplet.drawLine 180, 300, 180, 310
chartApplet.drawLine 230, 300, 230, 310
chartApplet.drawLine 280, 300, 280, 310
chartApplet.drawString "0          0.5          1          1.5          2          2.5          3", 28, 322
chartApplet.drawLine 30, 85,40, 85
chartApplet.drawLine 30, 160,40, 160
chartApplet.drawLine 30, 235,40, 235
chartApplet.drawString "0.02", 6, 15
chartApplet.drawString "0.01", 6, 164
chartApplet.drawString "0", 15, 315

'----
NM=2
NM = NM - 1
reDim Y(NM), Y0(NM), YW(NM), YS(NM), YQ(NM), YD(NM)
For  I=0 to 30
GZ(I)=0
GY(I)=0
GYI(I)=0
GX(I)=0
GXI(I)=0
next
'---- [a,b]
YA=0.0
YB=ValYB.value
'----
y1=0.02
yp2=Valyp2.value
y2=yp2/100
CL=ValCL.value
Y0(1)=0.00035
KYA = 309
KXA =5468
G = 47.0
L =1161*CL
ValL.value=L
x1=G*(y1-y2)/L
Valx1.value=x1
M = 26
D = -KXA / KYA
'----
Y0(0)=y1
GY(0)=y1
Y0(1)=x1
GZ(0)=0
GX(0)=x1
GXI(0)=(D*GX(0)-GY(0))/(D-M)
GYI(0)=M*GXI(0)
YH0=ValYH0.value
YE0=ValYE0.value
txtsolution.value="  充てん塔高さz   気相濃度 y   液相濃度 x" & vbCrLf
YH1 = YH0
X0 = YA

NZ=0
Do While Y0(0)>y2
NZ=NZ+1
FOR NK = 0 TO NM
YS(NK) = Y0(NK)
NEXT
YH = YH1
X = X0
Call  SUB1640
FOR NK = 0 TO NM
YQ(NK) = YW(NK)
NEXT
'1240------
YH = YH1 / 2.0
X = X0
Call SUB1640
FOR NK = 0 TO NM
YS(NK) = YW(NK)
NEXT
X = X0 + YH
Call SUB1640

FOR NK = 0 TO NM
	YD(NK) = (YW(NK) - YQ(NK))/15
	E1=Abs(YD(NK))/YH1/YW(NK)
	If E1>YE0 Then
		Valm.value="刻み幅が大きすぎ、終了"
		Exit Do
	End If
	If E1<(YE0/1000000) Then
		Valm.value="刻み幅が小さすぎ、終了"
		Exit Do
	End If
NEXT

'1500-----
X0 = X0 + YH1
FOR NK = 0 TO NM
Y0(NK) = YW(NK) + YD(NK)
Next
txtsolution.value=txtsolution.value & X0 &"    "& Y0(0) & "    " & Y0(1) & "    " & vbCrLf
GZ(NZ)=X0
GY(NZ)=Y0(0)
GX(NZ)=Y0(1)
GXI(NZ)=(D*GX(NZ)-GY(NZ))/(D-M)
GYI(NZ)=M*GXI(NZ)
Loop
Z=GZ(NZ-1)+(GZ(NZ)-GZ(NZ-1))*(GY(NZ-1)-y2)/(GY(NZ-1)-GY(NZ))
Valz.value=Z

chartApplet.setColor 255,0,0
For I=1 to NZ
chartApplet.drawLine GZ(I-1)*100+30, 310-GY(I-1)*15000, GZ(I)*100+30, 310-GY(I)*15000
NEXT
chartApplet.setColor 0,0,255
For I=1 to NZ
chartApplet.drawLine GZ(I-1)*100+30, 310-GX(I-1)*15000, GZ(I)*100+30, 310-GX(I)*15000
NEXT
chartApplet.setColor 0,255,0
For I=1 to NZ
chartApplet.drawLine GZ(I-1)*100+30, 310-GYI(I-1)*15000, GZ(I)*100+30, 310-GYI(I)*15000
NEXT
chartApplet.setColor 100,100,0
For I=1 to NZ
chartApplet.drawLine GZ(I-1)*100+30, 310-GXI(I-1)*15000, GZ(I)*100+30, 310-GXI(I)*15000
NEXT

END Sub

Sub Sub1640 
'------Runge-Kutta
FOR NK = 0 TO NM
Y(NK) = YS(NK)
YW(NK) = YS(NK)
NEXT
FOR NK = 0 TO NM
Call SUB2000
YDD = DYDX * YH
YW(NK) = YW(NK) + YDD / 6.0
Y(NK) = YS(NK) + YDD / 2.0
NEXT
X = X + YH / 2.0
FOR NK = 0 TO NM
Call SUB2000
YDD = DYDX * YH
YW(NK) = YW(NK) + YDD / 3.0
Y(NK) = YS(NK) + YDD / 2.0
NEXT
FOR NK = 0 TO NM
Call SUB2000
YDD = DYDX * YH
YW(NK) = YW(NK) + YDD / 3.0
Y(NK) = YS(NK) + YDD
NEXT
X = X + YH/2.0
FOR NK = 0 TO NM
Call SUB2000
YDD = DYDX * YH
YW(NK) = YW(NK) + YDD / 6.0
NEXT
End Sub

Sub SUB2000
'---------ﾋﾞﾌﾞﾝ ﾎｳﾃｲｼｷ
Select case NK
Case 0
' [ DYDX=Y(0)' ]
DYDX = -(KYA / G) * (Y(0) - M * (D * Y(1) - Y(0)) / (D - M))
Case 1
' [ DYDX=Y(1)' ]
DYDX = -(KXA / L) * ((D * Y(1) - Y(0)) / (D - M) - Y(1))
End Select
End Sub

--></SCRIPT>