計(jì)算機(jī)在化工中的應(yīng)用實(shí)驗(yàn)報(bào)告

1、計(jì)算機(jī)在化學(xué)化工中的應(yīng)用實(shí)驗(yàn)報(bào)告學(xué)院： 化學(xué)與化工學(xué)院 班級(jí)：12級(jí)碩勛勵(lì)志班 姓名: 徐凱杰 學(xué)號(hào)： 120702028實(shí)驗(yàn)一 傳熱實(shí)驗(yàn)中多變量的曲線的擬合一、實(shí)驗(yàn)?zāi)康?） 熟悉VB編程平臺(tái)2） 掌握多變量曲線擬合的算法3） 編擬合所給的傳熱實(shí)驗(yàn)?zāi)Ｐ偷腣B程序4） 通過(guò)實(shí)驗(yàn)數(shù)據(jù)求出模型數(shù)據(jù)、并掌握解線性方程組的克拉默法則二、運(yùn)行環(huán)境1） Microsoft Windows XP2） VB6.0三、實(shí)驗(yàn)原理略四、vb代碼Private Sub Command1_Click()Dim m As Integer'm=inputbox(“實(shí)驗(yàn)次數(shù)”)m = 7Dim x10, x20, y0

2、Dim i, j, k As IntegerDim a(1 To 10, 1 To 10), y(1 To 10), y1(1 To 10), a0, a1, a2Dim s, S1, S2, S3, b(1 To 10, 1 To 10), xxDim x1(1 To 10), x2(1 To 10), YY, sd'open"dem.dat"for input as#1'for i=1 to m' input#1,xx,YY' x1(i)=xx' x2(i0=xx2' y(i)=YY'next i'clos

3、e#1'7組努塞爾準(zhǔn)數(shù)、雷諾數(shù)及普蘭德準(zhǔn)數(shù)，數(shù)據(jù)最大時(shí)應(yīng)采用直接從文件讀取方法x10 = Array(0, 100, 200, 300, 500, 100, 700, 800) '注意下標(biāo)的起點(diǎn)處理（加0)x20 = Array(0, 2, 4, 1, 0.3, 5, 3, 4) '注意下標(biāo)的起點(diǎn)處理（加0)y0 = Array(0, 1.127, 2.416, 2.205, 2.312, 1.484, 6.038, 7.325) '注意下標(biāo)的起點(diǎn)處理（加0)For i = 1 To m x1(i) = Log(x10(i) x2(i) = Log(x20(i)

4、 y(i) = Log(y0(i)Next i'求解法方程系數(shù)矩陣a(1, 1) = ma(1, 2) = 0For i = 1 To m a(1, 2) = a(1, 2) + x1(i)Next ia(2, 1) = a(1, 2)a(1, 3) = 0For i = 1 To m a(1, 3) = a(1, 3) + x2(i)Next ia(3, 1) = a(1, 3)a(2, 2) = 0For i = 1 To m a(2, 2) = a(2, 2) + x1(i) * x1(i)Next ia(3, 3) = 0For i = 1 To m a(3, 3) = a(3

5、, 3) + x2(i) * x2(i)Next ia(2, 3) = 0For i = 1 To m a(2, 3) = a(2, 3) + x1(i) * x2(i)Next ia(3, 2) = a(2, 3)'求解法方程常數(shù)向量y1(1) = 0 For i = 1 To m y1(1) = y1(1) + y(i) Next iy1(2) = 0 For i = 1 To m y1(2) = y1(2) + x1(i) * y(i) Next iy1(3) = 0 For i = 1 To m y1(3) = y1(3) + x2(i) * y(i) Next i'(

6、利用克拉默法則解法方程/線性非常組)s = a(1, 1) * a(2, 2) * a(3, 3) + a(1, 2) * a(2, 3) * a(3, 1) + a(1, 3) * a(2, 1) * a(3, 2)s = s - a(1, 1) * a(2, 3) * a(3, 2) - a(1, 2) * a(2, 1) * a(3, 3) - a(1, 3) * a(2, 2) * a(3, 1)For j = 1 To 3 b(j, 1) = a(j, 1) a(j, 1) = y1(j)Next jS1 = a(1, 1) * a(2, 2) * a(3, 3) + a(1, 2)

7、 * a(2, 3) * a(3, 1) + a(1, 3) * a(2, 1) * a(3, 2)S1 = S1 - a(1, 1) * a(2, 3) * a(3, 2) - a(1, 2) * a(2, 1) * a(3, 3) - a(1, 3) * a(2, 2) * a(3, 1)For j = 1 To 3 a(j, 1) = b(j, 1)Next jFor j = 1 To 3 b(j, 2) = a(j, 2) a(j, 2) = y1(j) Next jS2 = a(1, 1) * a(2, 2) * a(3, 3) + a(1, 2) * a(2, 3) * a(3,

8、1) + a(1, 3) * a(2, 1) * a(3, 2)S2 = S2 - a(1, 1) * a(2, 3) * a(3, 2) - a(1, 2) * a(2, 1) * a(3, 3) - a(1, 3) * a(2, 2) * a(3, 1)For j = 1 To 3 a(j, 2) = b(j, 2) Next jFor j = 1 To 3 b(j, 3) = a(j, 3) a(j, 3) = y1(j)Next jS3 = a(1, 1) * a(2, 2) * a(3, 3) + a(1, 2) * a(2, 3) * a(3, 1) + a(1, 3) * a(2

9、, 1) * a(3, 2)S3 = S3 - a(1, 1) * a(2, 3) * a(3, 2) - a(1, 2) * a(2, 1) * a(3, 3) - a(1, 3) * a(2, 2) * a(3, 1)a0 = S1 / sa1 = S2 / sa2 = S3 / sText1.Text = Int(1000 * Exp(a0) + 0.5) / 1000 '四舍五入保留三位Text2.Text = Int(1000 * a1 + 0.5) / 1000Text3.Text = Int(1000 * a2 + 0.5) / 1000sd = 0For i = 1 T

10、o m sd = sd + Abs(a0 + a1 * x1(i) + a2 * x2(i) - y(i) '求 Nextsd = sd / mText4.Text = sd 'Int(1000 * sd + 0.5) / 1000Print Tab(50); "序號(hào)", "模型計(jì)算值", "實(shí)驗(yàn)值"For i = 1 To mPrintPrint Tab(45); i; (Text1.Text) * (x10(i) (Text2.Text) * (x20(i) (Text3.Text); 0.023 * (x10(i

11、) 0.8) * (x20(i) 0.3)NextEnd Sub五、實(shí)驗(yàn)結(jié)果截圖六、實(shí)驗(yàn)后思考。 VB編程是一種簡(jiǎn)單，并且效率高的可視化的、面向?qū)ο蠛筒捎檬录?qū)動(dòng)方式的結(jié)構(gòu)化高級(jí)程序設(shè)計(jì)語(yǔ)言。通過(guò)對(duì)本實(shí)驗(yàn)的實(shí)際操作，我掌握了多變量曲線擬合的基本算法，了解了解線性方程組的克拉默法則。并且，同時(shí)在以后的工作中，可以通過(guò)這個(gè)實(shí)驗(yàn)來(lái)解決大部分實(shí)驗(yàn)數(shù)據(jù)及模型參數(shù)的擬合問(wèn)題。實(shí)驗(yàn)二 梯度法擬合蒸汽壓與溫度關(guān)系模型一 、實(shí)驗(yàn)?zāi)康?) 掌握梯度法擬合的基本算法以及理解其普適性2) 編寫梯度法擬合蒸汽壓與溫度的關(guān)系的VB程序3) 通過(guò)實(shí)對(duì)程序進(jìn)行驗(yàn)證，并注意比較初值對(duì)運(yùn)行速度和結(jié)果的影響二 、運(yùn)行環(huán)境1） Mi

12、crosoft Windows XP2） VB6.0三 、實(shí)驗(yàn)原理略四、 實(shí)驗(yàn)VB程序代碼Private Sub Command1_Click(Index As Integer)Dim m, n As Integerm = 6Dim i, j, k As IntegerDim A, B, C, F, ee, P(1 To 10), T(1 To 10)Dim A1, B1, C1, TA, TB, TC, TT, f1, f2, f3Dim sd, W, S, EY, XX, YY'(由dem.dat輸入實(shí)驗(yàn)數(shù)據(jù)XX = Array(-23.7, -10, 0, 10, 20, 30,

13、 40) '注意下標(biāo)的起點(diǎn)處理（加0)YY = Array(0.101, 0.174, 0.254, 0.359, 0.495, 0.662, 0.88) '注意下標(biāo)的起點(diǎn)處理（加0)Print "直接讀數(shù)據(jù)文件后計(jì)算"For i = 1 To m T(i) = XX(i) T(i) = 273.15 + T(i) P(i) = YY(i) * 7600 Print T(i), P(i)Next iClose i A = Val(InputBox("A") '指定初值 B = Val(InputBox("B")

14、 '指定初值 C = Val(InputBox("C") '指定初值1000 F = 0For i = 1 To m ee = FNP(A, B, C, T(i), P(i) ee = ee 2 F = F + eeNext if1 = 0A1 = A + 0.000001 * A'print"A,A1="A,A1For i = 1 To m ee = FNP(A1, B, C, T(i), P(i) ee = ee 2 f1 = f1 + eeNext iTA = (f1 - F) / (0.000001 * A)'pr

15、int f1,F,TA'A=val(inputbox("A")f2 = 0B1 = B + 0.00001 * BFor i = 1 To m ee = FNP(A, B1, C, T(i), P(i) ee = ee 2 f2 = f2 + eeNext iTB = (f2 - F) / (0.00001 * B)f3 = 0C1 = C + 0.00001 * CFor i = 1 To m ee = FNP(A, B, C1, T(i), P(i) ee = ee 2 f3 = f3 + eeNext iTC = (f3 - F) / (0.00001 * C

16、)TT = TA 2 + TB 2 + TC 2TT = Sqr(TT)If TT > 0.001 ThenA = A - 0.005 * TAB = B - 1.5 * TBC = C - 0.001 * TCGoTo 1000ElseEnd IfPrintsd = 0For i = 1 To m '/計(jì)算絕對(duì)平均相對(duì)誤差sd = sd + Abs(FNSD(A, B, C, T(i), P(i) / P(i)Print FNSD(A, B, C, T(i), P(i)Next isd = sd / mPrintPrint "A,B,C=" A, B, CP

17、rint "sd=" sd '/打印絕對(duì)平均相對(duì)誤差End SubPublic Function FNP(A, B, C, T, P)FNP = (A - B / (T + C) - Log(P)End FunctionPublic Function FNSD(A, B, C, T, P)FNSD = Exp(A - B / (T + C) - PEnd Function五 、實(shí)驗(yàn)結(jié)果截圖六 、實(shí)驗(yàn)后思考。本實(shí)驗(yàn)是基于最小二乘原理，函數(shù)擬合的目標(biāo)是使擬合函數(shù)和實(shí)際測(cè)量值之間的差的平方和最小。對(duì)于最小值的問(wèn)題，梯度法是用負(fù)梯度方向作為優(yōu)化搜索方向。而梯度法是一個(gè)簡(jiǎn)單的

18、迭代優(yōu)化計(jì)算方法。注意的是，負(fù)梯度的最速下降性是一個(gè)局部的性質(zhì)。在計(jì)算的前期使用此法，當(dāng)接近極小點(diǎn)時(shí)，在改用其他的算法，如共軛梯度法。 實(shí)驗(yàn)三 二分法求解化工中的非線性方程一、實(shí)驗(yàn)?zāi)康?） 掌握二分法解非線性方程組的基本算法2） 編寫二分法邱玨非線性方程組的VB程序3） 通過(guò)實(shí)例的程序進(jìn)行調(diào)試，并學(xué)習(xí)輸出數(shù)據(jù)格式化二、 運(yùn)行環(huán)境1） Microsoft WindowsXP2） VB6.0三 、實(shí)驗(yàn)原理略四 、實(shí)驗(yàn)VB代碼Private Sub Command1_Click()Dim ax As SingleDim bx As SingleDim cx As SingleDim ay As Si

19、ngleDim by As SingleDim cy As SingleDim e As SingleDim num As Integer '累計(jì)次數(shù)變量Dim st As StringDim ch As StringDim sp As Stringch = Chr(13) + Chr(10)sp = Space(10)st = "二分法解方程" + chst = st + "求2,3-二甲基苯胺沸點(diǎn)（當(dāng) P=101325 時(shí) 解 lnP=59.7622-8013.69/T-5.081lnT)" + chax = 200bx = 500e = 0

20、.01st = st + "區(qū)間左端點(diǎn)初始值 ax=" + Str(ax) + chst = st + "區(qū)間右端點(diǎn)初始值 bx=" + Str(bx) + chst = st + "精度控制限 e=" + Str(e) + chst = st + "num" + sp + "ax" + Space(14) + "bx" + Space(14) + "|ax-bx|" + chay = F(ax)by = F(bx)num = 1Do While Abs(

21、ax - bx) > e cx = (ax + bx) / 2 cy = F(cx) If cy = 0 Then Exit Do '如果已得解，則退出循環(huán) If cy * ay > 0 Then ax = cx ay = cy Else bx = cx by = cy End If st = st + Format(num, "000") + sp + Format(ax, "000.00") + sp + Format(bx, "000.00") + sp + Format(Abs(ax - bx), &quo

22、t;0.0000") + ch num = num + 1Loopst = st + ch + "2,3-二甲基苯胺沸點(diǎn):" + Format(cx, "#00.00") + "K" + ch + chst = st + "*時(shí)間:" + Str(Time) + Space(3) + "日期:" + Str(Date) + chText1.Text = ""Text1.Text = stEnd Sub'二分法求2,3-二甲基苯胺沸點(diǎn)所用函數(shù)Private F

23、unction F(ByVal u As Single)F = Log(101325) - 59.7622 + 8013.69 / u + 5.081 * Log(u) '注意對(duì)數(shù)運(yùn)算End Function五、實(shí)驗(yàn)結(jié)果截圖六 、實(shí)驗(yàn)后思考。 通過(guò)應(yīng)用微積分中的介值定理，是是用二分法的前提條件。如果我們所要求解的方程從物理意義上來(lái)講確實(shí)存在實(shí)根，但又不滿足f(a)f(b)<0,這時(shí)候，我們必須通過(guò)改變a和b的值來(lái)滿足二分法的應(yīng)用條件。實(shí)驗(yàn)四 主元最大高斯消元法解化工中的線性方程組一、實(shí)驗(yàn)?zāi)康?） 掌握主元最大高斯消元法2） 編寫最大高斯消元法求解線性方程組的VB程序3） 通過(guò)實(shí)例

24、對(duì)程序進(jìn)行調(diào)試，并比較一般的高斯消去法比較二 、運(yùn)行環(huán)境1） Microsoft WIndowsXP_2） VB6.0三 、實(shí)驗(yàn)原理略四 、實(shí)驗(yàn)程序代碼Private Sub Command1_Click()Dim m, n As IntegerDim a(), z(), x(), w, aa(), s, t, k, ln = 4ReDim a(n + 2, 2 + n), z(n + 2, 2 + n), x(n + 1), aa(n + 2, 2 + n)Dim i, j, k1, k2, stDim ch As StringDim sp As Stringch = Chr(13) + C

25、hr(10)sp = Space(5)a(1, 1) = 6# / 123.1a(1, 2) = 6# / 93.13a(1, 3) = 3# / 73.1a(1, 4) = 2# / 43.07a(2, 1) = 5# / 123.1a(2, 2) = 7# / 93.13a(2, 3) = 7# / 73.1a(2, 4) = 6# / 43.07a(3, 1) = 1# / 123.1a(3, 2) = 1# / 93.13a(3, 3) = 1# / 73.1a(3, 4) = 0# / 43.07a(4, 1) = 2# / 123.1a(4, 2) = 0# / 93.13a(4,

26、 3) = 1# / 73.1a(4, 4) = 1# / 43.07a(1, 5) = 57.78 / 12.01a(2, 5) = 7.92 / 1.008a(3, 5) = 11.23 / 14.01a(4, 5) = 23.09 / 16st = st + "主元最大高斯消去法解線性方程組" + chst = st + "設(shè)有一混合物由硝基苯、苯胺、氨基丙酮、乙醇組成;" + chst = st + "對(duì)該混合物進(jìn)行元素分析結(jié)果以百分?jǐn)?shù)表示如下" + chst = st + "C%=57.78%;H%=7.92%;N

27、%=11.23%;O%=23.09%" + chst = st + "原子量：A(C)=12.01;A(H)=1.008;A(N)=14.01;A(O)=16.00" + chst = st + "分子量：硝基苯 123.1;苯胺 93.13;氨基丙酮 73.10;乙醇 43.07" + chst = st + "硝基苯分子C-6;H-5;N-1;O-2" + chst = st + "苯胺分子C-6;H-7;N-1;O-0" + chst = st + "氨基丙酮分子C-3;H-7;N-1;O

28、-1" + chst = st + "乙醇分子C-2;H-6;N-0;O-1" + chst = st + "確定上面四種化合物在混合物中所占的百分比" + ch + ch'尋找主元For i = 1 To n If i = n Then GoTo 200 For t = i + 1 To n If Abs(a(i, i) < Abs(a(t, i) Then For s = i To n + 1 aa(t, s) = a(i, s) a(i, s) = a(t, s) a(t, s) = aa(t, s) Next s Else

29、 End If Next t200'消去w = a(i, i) For j = 1 To n + 1 a(i, j) = a(i, j) / w Next jIf i = n Then GoTo 100For j = i + 1 To n For k = i + 1 To n + 1 z(i, k) = a(i, k) * a(j, i) a(j, k) = a(j, k) - z(i, k) Next kNext jNext i100'回代x(n + 1) = 0 For k = n To 1 Step -1 s = 0 For j = k + 1 To n s = s +

30、a(k, j) * x(j) Next j x(k) = a(k, n + 1) - s 'st=st+"x("+str(i)+")="+format(x(i),"00.00")+"%"+ch 'print"x("k;")="x(k) Next kFor i = 1 To n '輸出結(jié)果 st = st + "x(" + Str(i) + ")=" + Format(x(i), "00.00"

31、) + "%" + chNext ist = st + chst = st + "*時(shí)間:" + Str(Time) + Space(3) + "日期:" + Str(Date) + chText1.Text = ""Text1.Text = stEnd Sub五 、實(shí)驗(yàn)結(jié)果截圖六、 實(shí)驗(yàn)后思考高斯消去法不需要方程組的初值，也不需要重復(fù)迭代計(jì)算。只通過(guò)“消去”和“回代”2個(gè)過(guò)程就可以直接求出方程組的解。然后若是在消去的過(guò)程中，若碰到主元為0，則無(wú)法計(jì)算。所以，發(fā)展了“主元最大高斯消去法”。就是在主元所在的列中，尋找

32、到最大的元素，進(jìn)行行與行之間的調(diào)換，并將該最大的元素作為主元，保證主元不為0。實(shí)驗(yàn)五 松弛迭代法求解化工中的線性方程組一、實(shí)驗(yàn)?zāi)康?） 掌握松弛迭代法的基本算法及和緊湊迭代的細(xì)微區(qū)別2） 編寫松弛迭代法求救線性方程組的VB代碼，注意學(xué)習(xí)從文件讀取數(shù)據(jù)3） 通過(guò)實(shí)例的程序進(jìn)行驗(yàn)證，并觀察松弛迭代因子對(duì)結(jié)果的影響二、運(yùn)行環(huán)境1） Microsoft WIndowse XP2） VB6.0三、實(shí)驗(yàn)原理略四、實(shí)驗(yàn)程序代碼Private Sub Command1_Click()Dim n As IntegerDim i, j, ff, t, k, l, hDim st As StringDim a()

33、As SingleDim y() As SingleDim b() As SingleDim g() As SingleDim x1() As SingleDim x2() As SingleDim jk() As IntegerDim ch As StringDim sp As Stringch = Chr(13) + Chr(10)sp = Space(5)CommonDialog1.CancelError = True'on error goto errhandlerCommonDialog1.Filter = "數(shù)據(jù)文件（*.txt）|*.txt|拉圖文件（*.bmp

34、）|*.bmp|AllFiles(*.*)|*.*" '文件過(guò)濾CommonDialog1.FilterIndex = 0CommonDialog1.DialogTitle = "加載增廣矩陣數(shù)據(jù)文件"CommonDialog1.ShowOpen'*數(shù)據(jù)文件的行數(shù)就是方程的個(gè)數(shù)Open CommonDialog1.FileName For Input As #1Do While Not EOF(1) Line Input #1, st n = n + 1LoopClose #1'*數(shù)據(jù)文件的行數(shù)就是方程的個(gè)數(shù)ReDim a(1 To n,

35、1 To n) As SingleReDim b(1 To n, 1 To n) As SingleReDim x1(1 To n) As SingleReDim x2(1 To n) As SingleReDim g(1 To n) As SingleReDim y(1 To n) As Single'*讀數(shù)據(jù)Open CommonDialog1.FileName For Input As #1For i = 1 To n For j = 1 To n Input #1, a(i, j) '方程等號(hào)左端數(shù)據(jù) Next j Input #1, y(i) '方程等號(hào)右端數(shù)

36、據(jù)Next iClose #q'*讀數(shù)據(jù)st = "松弛迭代法解線性方程" + chst = st + Space(5) + "數(shù)據(jù)來(lái)源于" + CommonDialog1.FileName + ch + chst = st + "增廣矩陣如下（對(duì)二甲苯-間二甲苯-鄰二甲苯-乙苯-（混合物）：" + ch + chst = st + Space(5) + "第一行為12.5nm波長(zhǎng)處摩爾吸收系數(shù)-混合物吸收" + chst = st + Space(5) + "第一行為13.0nm波長(zhǎng)處摩爾吸收系

37、數(shù)-混合物吸收" + chst = st + Space(5) + "第一行為13.4nm波長(zhǎng)處摩爾吸收系數(shù)-混合物吸收" + chst = st + Space(5) + "第一行為14.3nm波長(zhǎng)處摩爾吸收系數(shù)-混合物吸收" + ch + ch'*輸出原始數(shù)據(jù) For i = 1 To n For j = 1 To n If a(i, j) >= 0 Then st = st + Space(5) + Format(a(i, j), "0.00000") Else st = st + Space(4) +

38、Format(a(i, j), "0.00000") End If Next j If y(i) >= 0 Then st = st + Space(5) + Format(y(i), "0.00000") + ch Else st = st + Space(4) + Format(y(i), "0.00000") + ch End If Next i'*輸出原始數(shù)據(jù)'-For i = 1 To n x1(i) = 0 x2(i) = 0Next i'for i = 1 to n' for j =

39、 1 to n' a(i,j) = InputBox("a("&i&","&j&")")' Print a(i,j),' next j'y(i) = InputBox("y("&i&")")'print" ",y(i)'Next i'產(chǎn)生迭代矩陣For i = 1 To n g(i) = y(i) / a(i, i) For j = 1 To n If j = i Then

40、 b(i, j) = 0 Else b(i, j) = -a(i, j) / a(i, i) End If Next j Next i e = InputBox("輸入松弛因子") '開始松弛迭代Do If k >= 1 Then For i = 1 To n x1(i) = x2(i) Next i End If For i = 1 To n s = g(i) For j = 1 To n s = s + b(i, j) * x2(j) Next j x2(i) = (1 - e) * x1(i) + e * s '注意 Next i eer = 0

41、 For i = 1 To n eer = cer + Abs(x1(i) - x2(i) '計(jì)算誤差 Next i k = k + 1 '累計(jì)次數(shù)Loop While (k < 100 And eer >= 0.001) Print kst = st + ch + "方程組的解為：" + ch + chFor i = 1 To n st = st + "x(" + Str(i) + ")=" + Format(x2(i), "0.00000") + ch Next ist = st +

42、 ch + "迭代次數(shù)為：" + Str(k) + ch 'format(k,"000")st = st + ch + "松弛因子為：" + Format(e, "0.0000") + chst = st + ch + "誤差為：" + Format(eer, "0.000000") + chst = st + chst = st + "*時(shí)間：" + Str(Time) + Space(3) + "日期：" + Str(Dat

43、e) + chText1.Text = ""Text1.Text = st End Sub五、實(shí)驗(yàn)結(jié)果截圖六、 實(shí)驗(yàn)后思考。松弛迭代法是數(shù)值計(jì)算中解線性代數(shù)方程組的一類迭代法。逐次超松弛迭代過(guò)程中，已知迭代方程及其系數(shù)矩陣，對(duì)任意的初始值，確定超松弛因子，用迭代矩陣來(lái)進(jìn)行計(jì)算確定譜半徑，然后其絕對(duì)值小于一解出來(lái)超松弛因子。而緊湊迭代是當(dāng)松弛因子為1的時(shí)候，叫做緊湊迭代。兩者的區(qū)別在于松弛因子的不同。實(shí)驗(yàn)六 龍格庫(kù)塔法求解化工過(guò)程中的常微分方程一、實(shí)驗(yàn)?zāi)康?） 掌握龍格庫(kù)塔法的基本原理2） 編寫龍格庫(kù)塔法解決常微分方程的VB程序3） 通過(guò)實(shí)例的程序進(jìn)行調(diào)試和驗(yàn)證，并觀察初值對(duì)

44、計(jì)算過(guò)程及結(jié)果的影響4） 掌握VB繪制二維曲線圖的方法和繪圖參數(shù)的設(shè)置二、 運(yùn)行環(huán)境1） Microsoft WindowsXP2） VB6.0三 、實(shí)驗(yàn)原理略四、 實(shí)驗(yàn)程序截圖Private Sub Command1_Click()Const eps = 0.00001Dim t() As SingleDim x() As SingleDim y() As SingleDim z() As SingleDim J1, J2 As SingleDim K1, K2, K3, K4 As SingleDim Q1, Q2, Q3, Q4 As SingleDim S1, S2, S3, S4 A

45、s SingleDim h As SingleDim i As IntegerDim n As Integerh = 0.01J1 = 1J2 = 1.1n = Int(10 / h)ReDim t(n + 1), x(n + 1), y(n + 1), z(n + 1) As Singlet(0) = 0x(0) = 0y(0) = 0'z(0) = 0For i = 0 To n - 1K1 = -J1 * x(i)Q1 = J1 * x(i) - J2 * y(i)S1 = h * (J2 * y(i)K2 = -J1 * (x(i) + h * K1 / 2)Q2 = J1 *

46、 (x(i) + h * K1 / 2) - J2 * (y(i) + h * Q1 / 2)S2 = h * (J2 * (y(i) + Q1 / 2)K3 = -J1 * (x(i) + h * K2 / 2)Q3 = J1 * (x(i) + h * K2 / 2) - J2 * (y(i) + h * Q2 / 2)S3 = h * (J2 * (y(i) + Q2 / 2)K4 = -J1 * (x(i) + h * K3)Q4 = J1 * (x(i) + h * K3) - J2 * (y(i) + h * Q3)S4 = h * (J2 * (y(i) + Q3)x(i + 1

47、) = x(i) + h * (K1 + 2 * K2 + 2 * K3 + K4) / 6 '計(jì)算A物質(zhì)的濃度y(i + 1) = y(i) + h * (Q1 + 2 * Q2 + 2 * Q3 + Q4) / 6z(i + 1) = z(i) + (S1 + 2 * Q2 + 2 * Q3 + Q4) / 6z(i) = x(0) - x(i) - y(i)t(i + 1) = t(i) + h '計(jì)算反應(yīng)時(shí)間tNext iDim axisname1 As String, axisname2 As Stringaxisname1 = "t"axisna

48、me2 = ""xy_axis picture1, t(), x(), axisname1, axisname2For i = 0 To n - 1 picture1.PSet (t(i), x(i)Next i'xy_axis picture1,t(),y(),axisname1,axisname2For i = 0 To n - 1 Picture.PSet (t(i), y(i) Next i 'xy_axis picture1,t(),z(),axisname1,axisname2 For i = 0 To n - 1 Picture.PSet (t

49、(i), z(i)Next iEnd Sub模塊代碼Sub xy_axis(pic As PictureBox, x1() As Single, y1() As Single, axisname1 As String, axisname2 As String)On Error GoTo problemx:Dim maxnumber As Single, minnumber As SingleDim leftx As Single, topy As SingleDim rightx As Single, bottomy As SingleDim n As Integern = UBound(x1)pic.Font