Додаток І. Програма розрахунку параметрів рівноваги газофазних оборотних реакцій

І.1. Програма на мові BASIC

10 REM  ********************************************************

20 REM  *              Программа Р01 RAVN                      *

30 REM  ********************************************************

40 REM  *Расчет параметров состояния равновесия широкого класса*

50 REM  *            газофазных обратимых реакций              *

60 REM  *                                                      *

70 REM  *               aA + bB <=====> rR + sS                *

80 REM  *                                                      *

90 REM  ********************************************************

100 REM *          С ИСПОЛЬЗОВАНИЕМ МЕТОДА НЬЮТОНА             *

110 REM ********************************************************

120 OPTION BASE 1

130 DIM N(4), A(4), FU(8), B(4), P(4), C(4)

140 PRINT "****************************************************"

150 PRINT "*                                                  *"

160 PRINT "*                РЕШЕНИЕ ЗАДАЧИ                    *"

170 PRINT "*                                                  *"

180 PRINT "****************************************************"

190 PRINT

200 PRINT "************ВВОД ИСХОДНЫХ ДАННЫХ**************": PRINT

205 PRINT "ВВЕДИТЕ:": PRINT

210 INPUT "NA, NB, NR, NS - Стехиoметрические коэффициенты "; N(1), N(2), N(3), N(4)

220 IF N(1) + N(2) = 0 THEN STOP

230 PRINT "A0, B0, R0, S0 - исходные количества"

235 INPUT "                 компонентов, Кмоль             "; A(1), A(2), A(3), A(4)

240 INPUT "KP - константа равновесия, МПа^(r+s-a-b)        "; RP

250 INPUT "T - температура, К                              "; T

260 INPUT "P - давление, МПа                               "; P

270 INPUT "X - нач. приближение характеристик равновесия   "; X

280 REM

290 REM Определение констант

300 REM ITMAX - Допустимое число операций

310 REM CV, DV - критерии сходимости

320 REM

330 ITMAX = 20

340 CV = .0001

350 DV = 100000!

360 NDELT = N(3) + N(4) - (N(1) + N(2)):

370 NSM0 = 0: FOR I = 1 TO 4

380 NSM0 = NSM0 + A(I): NEXT I: VM = 22.4 * T * .1 / 273 / P: V0 = NSM0 * VM

390 KP = RP / (8.310001E-03 * T) ^ NDELT * VM ^ NDELT

400 PRINT

410 PRINT "****************РЕЗУЛЬТАТЫ РАСЧЕТА*************": PRINT

420 PRINT SPC(10); "Последовательные приближения"

430 PRINT

450 PRINT SPC(2); "I"; TAB(12); "X"; TAB(28); "DX"; TAB(44); "F(X)"

460 PRINT

470 REM

480 REM Итерационный цикл

490 REM

500 FOR II = 1 TO ITMAX

510 FOR J = 1 TO 4

520 IPR = N(J): IGO = J: GOSUB 1050: FU(J) = FUNC: NEXT J

530 FOR J = 5 TO 8

540 IPR = N(J - 4) - 1: IGO = J - 4: GOSUB 1050: FU(J) = FUNC: NEXT J

550 NSMR = NSM0 + NDELT * X

560 FO = KP * NSMR ^ NDELT * FU(1) * FU(2) - FU(3) * FU(4)

570 DF1 = -KP * NSMR ^ NDELT * (N(2) * N(2) * FU(1) * FU(6) + N(1) * N(1) * FU(2) * FU(5))

580 DF2 = KP * NDELT ^ 2 * NSMR ^ (NDELT - 1) * FU(1) * FU(2)

590 DF0 = DF1 + DF2 - N(4) * N(4) * FU(3) * FU(8) - N(3) * N(3) * FU(4) * FU(7)

600 IF ABS(DF0) <= 1E-25 GOTO 730

610 DX = FO / DF0

620 PRINT USING "##     #.#####^^^^"; II; X;

630 PRINT USING "   #.#####^^^^"; DX; FO

640 REL = ABS(DX / X)

650 IF REL <= CV GOTO 810

660 IF REL >= DV GOTO 750

670 X = X - DX: NEXT II

680 REM

690 REM ПЕЧАТЬ СООБЩЕНИЙ ОБ ОШИБКЕ И ВОЗВРАТ К НАЧАЛУ ПРОГРАММЫ

700 REM

710 PRINT " НЕТ СХОДИМОСТИ ПРИ ITMAX="; ITMAX

720 GOTO 210

730 PRINT " ФУНКЦИЯ DF(X) ПРИБЛИЖАЕТСЯ К 0"

740 GOTO 210

750 PRINT "МЕТОД РАСХОДИТСЯ"

760 GOTO 210

770 REM

780 REM РАСЧЕТ РАВНОВЕСНЫХ КОЛИЧЕСТВ, РАВНОВЕСНЫХ ПАРЦИАЛЬНЫХ

790 REM ДАВЛЕНИЙ КОМПОНЕНТОВ И КОНСТАНТЫ РАВНОВЕСИЯ

800 REM

810 FOR I = 1 TO 4

820 IF I <= 2 THEN B(I) = A(I) - N(I) * X ELSE B(I) = A(I) + N(I) * X

830 NEXT I

840 V = VM * NSMR

850 FOR I = 1 TO 4

860 P(I) = B(I) * 8.310001E-03 * T / V: C(I) = B(I) / NSMR: NEXT I

870 PSUM = 0!: FOR I = 1 TO 4: PSUM = PSUM + P(I): NEXT I

880 PRINT

890 PRINT SPC(10); "РАВНОВЕСНЫЕ КОЛИЧЕСТВА КОМПОНЕНТОВ": PRINT

900 PRINT "AEQ="; B(1); "BEQ="; B(2); "REQ="; B(3); "SEQ="; B(4); "[кмоль]"

910 PRINT : PRINT SPC(5); "ОБЪЕМ НАЧАЛЬНОЙ СМЕСИ V0="; V0; "[м^3]": PRINT

920 PRINT SPC(5); "ОБЪЕМ РАВНОВЕСНОЙ СМЕСИ V="; V; "[м^3]": PRINT

930 PRINT SPC(5); "РАВНОВЕСНЫЕ ПАРЦИАЛЬНЫЕ ДАВЛЕНИЯ КОМПОНЕНТОВ"

940 PRINT : PRINT "PA="; P(1); "PB="; P(2); "[МПа]": PRINT

950 PRINT SPC(5); "РАВНОВЕСНОЕ ОБЩЕЕ ДАВЛЕНИЕ P="; PSUM; "[МПа]": PRINT

960 PRINT SPC(5); "РАВНОВЕСНЫЕ МОЛЯРНЫЕ ДОЛИ КОМПОНЕНТОВ": PRINT

970 PRINT "CA="; C(1); "CB="; C(2); "CR="; C(3); "CS="; C(4): PRINT

980 FOR I = 1 TO 4

990 IF P(I) = 0! THEN P(I) = 1!

995 NEXT I

1000 CALCRP = P(3) ^ N(3) * P(4) ^ N(4) / (P(1) ^ N(1) * P(2) ^ N(2))

1010 PRINT "РАСЧИТАННАЯ КР="; CALCRP; "[МПа ^ ("; NDELT; ")]"

1020 PRINT : PRINT "ВЫПОЛНЕНИЕ ПРОГРАММЫ ЗАВЕРШЕНО"

1030 END

1040 REM

1050 REM ПОДПРОГРАММА FUNC

1070 REM

1080 IF IGO > 2 THEN K = 1 ELSE K = -1

1090 IF N(IGO) <= 0 GOTO 1110 ELSE BASE1 = A(IGO) + K * N(IGO) * X

1100 IF BASE1 > 0! THEN FUNC = BASE1 ^ IPR: RETURN

1110 FUNC = 1!

1120 RETURN

1130 END

 

 

 

 

 

І.2. Програма на мові FORTRAN

      REAL N,ITMAX,NDELT,NSM0,II,NSMR,I,IPR,IGO,J,K,KP

      DIMENSION N(4), A(4), FU(8), B(4), P(4), C(4)

      WRITE(*,*) '****************************************************'

      WRITE(*,*) '*                РЕШЕНИЕ ЗАДАЧИ                    *'

      WRITE(*,*) '****************************************************'

      WRITE(*,*) '************ВВОД ИСХОДНЫХ ДАННЫХ**************'

      WRITE(*,*)'ВВЕДИТЕ:'

5     WRITE(*,*)'NA, NB, NR, NS - Стехиoметрические коэффициенты'

      READ(*,*)N(1),N(2),N(3),N(4)

      ZZ=N(1)+N(2)

      IF(ZZ.EQ.0) STOP

      WRITE(*,*)'A0,B0,R0,S0 - исходные количества компонентов, Кмоль'

      READ(*,*)A(1),A(2),A(3),A(4)

      WRITE(*,*)'KP - константа равновесия, МПа^(r+s-a-b)'

      READ(*,*)RP

      WRITE(*,*)'T - температура, К'

      READ(*,*)T

      WRITE(*,*)'P - давление, МПа'

      READ(*,*)PE

      WRITE(*,*)'X - нач. приближение характеристик равновесия'

      READ(*,*)X

      ITMAX=40.

      CV=.0001

      DV=100000.

      NDELT = N(3)+N(4)-(N(1)+N(2))

      NSM0=0.

      DO I=1,4

      NSM0=NSM0+A(I)

      END DO

      VM=22.4*T*.1/(273.*PE)

      V0=NSM0*VM

      KP=RP/(((0.008310001*T)**NDELT)*(VM**NDELT))

      WRITE(*,*)

      WRITE(*,*)'****************РЕЗУЛЬТАТЫ РАСЧЕТА*************'

      WRITE(*,*)'          Последовательные приближения'

      WRITE(*,*)'  I          X        DX                        F(X)'

C     Итерационный цикл

      II=1

9     DO J=1,4

      IPR=N(J)

      IGO=J

      FU(J)=FUNC(IGO,K,IPR,X,a,n)

      END DO

      DO J=5,8

      IPR=N(J-4)-1

      IGO=J-4

      FU(J)=FUNC(IGO,K,IPR,X,a,n)

      END DO

      NSMR=NSM0+NDELT*X

      FO=KP*NSMR**NDELT*FU(1)*FU(2)-FU(3)*FU(4)

      DF1=-KP*NSMR**NDELT*(N(2)*N(2)*FU(1)*FU(6)+N(1)*N(1)*FU(2)*FU(5))

      DF2=KP*NDELT**2*NSMR**(NDELT-1)*FU(1)*FU(2)

      DF0=DF1+DF2-N(4)*N(4)*FU(3)*FU(8)-N(3)*N(3)*FU(4)*FU(7)

      IF (ABS(DF0).LE.1/(10**25))GOTO 30

      DX=FO/DF0

      WRITE(*,10)II,X

10    FORMAT(F7.3,F7.2)

      WRITE(*,100)DX,FO

100   format(f7.3,f20.0)

      REL=ABS(DX/X)

      IF(REL.LE.CV) GOTO 40

      IF(REL.GE.DV) THEN

      GOTO 35

      X=X-DX

      END IF

      IF(II.EQ.ITMAX) THEN

      GOTO 15

      ELSE

      II=II+1

      GOTO 9

15    END IF

C     ПЕЧАТЬ СООБЩЕНИЙ ОБ ОШИБКЕ И ВОЗВРАТ К НАЧАЛУ ПРОГРАММЫ

      WRITE(*,*)'НЕТ СХОДИМОСТИ ПРИ ITMAX=',ITMAX

      GOTO 5

30    WRITE(*,*)'ФУНКЦИЯ DF(X) ПРИБЛИЖАЕТСЯ К 0'

      GOTO 5

35    WRITE(*,*)'МЕТОД РАСХОДИТСЯ'

      GOTO 5

C     РАСЧЕТ РАВНОВЕСНЫХ КОЛИЧЕСТВ, РАВНОВЕСНЫХ ПАРЦИАЛЬНЫХ

C     ДАВЛЕНИЙ КОМПОНЕНТОВ И КОНСТАНТЫ РАВНОВЕСИЯ

40    CONTINUE

      DO I=1,4

      IF(I.LE.2) THEN

      B(I)=A(I)-N(I)*X

      ELSE

      B(I)=A(I)+N(I)*X

      END IF

      END DO

      V=VM*NSMR

      DO I=1,4

      P(I)=B(I)*.00831000188*T/V

      C(I)=B(I)/NSMR

      END DO

      PSUM=0

      DO I=1,4

      PSUM=PSUM+P(I)

      END DO

      WRITE(*,*)

      WRITE(*,*)'          РАВНОВЕСНЫЕ КОЛИЧЕСТВА КОМПОНЕНТОВ'

      WRITE(*,*)

      WRITE(*,*)'AEQ=',B(1),'BEQ=',B(2),'REQ=',B(3),'SEQ=',B(4),

     * '[кмоль]'

      WRITE(*,*)

      WRITE(*,*)'     ОБЪЕМ НАЧАЛЬНОЙ СМЕСИ V0=',V0,'[м^3]'

      WRITE(*,*)

      WRITE(*,*)'     ОБЪЕМ РАВНОВЕСНОЙ СМЕСИ V=',V,'[м^3]'

      WRITE(*,*)

      WRITE(*,*)'     РАВНОВЕСНЫЕ ПАРЦИАЛЬНЫЕ ДАВЛЕНИЯ КОМПОНЕНТОВ'

      WRITE(*,*)

      WRITE(*,*)'PA=',P(1),'PB=',P(2),'[МПа]'

      WRITE(*,*)

      WRITE(*,*)'     РАВНОВЕСНОЕ ОБЩЕЕ ДАВЛЕНИЕ P=',PSUM,'[МПа]'

      WRITE(*,*)

      WRITE(*,*)'     РАВНОВЕСНЫЕ МОЛЯРНЫЕ ДОЛИ КОМПОНЕНТОВ'

      WRITE(*,*)

      WRITE(*,*)'CA=',C(1),'CB=',C(2),'CR=',C(3),'CS=',C(4)

      WRITE(*,*)

      DO I=1,4

      IF(P(I).EQ.0) P(I)=1

      END DO

      CALCRP=P(3)**N(3)*P(4)**N(4)/(P(1)**N(1)*P(2)**N(2))

      WRITE(*,*)'РАСЧИТАННАЯ КР=',CALCRP,'[МПа^(',NDELT,')]'

      WRITE(*,*)

      WRITE(*,*)'ВЫПОЛНЕНИЕ ПРОГРАММЫ ЗАВЕРШЕНО'

      END

     

      FUNCTION FUNC(IGO,K,IPR,X,a,n)

      dimension a(4),n(4)

      REAL IGO,K,IPR,n

      IF (IGO.GT.2) THEN

      K=1

      ELSE

      K=-1

      END IF

      IF(N(IGO).LE.0) THEN

      GOTO 111

      ELSE

      BASE1=a(igo)+K*n(igo)*X

      END IF

      IF (BASE1.gt.0)    THEN

      FUNC=BASE1**IPR

      RETURN

      ELSE

111   FUNC=1

      RETURN

      END IF

      END

     

     

І.3. Програма на мові PASCAL

uses math;

label l1,l2;

const

  ITMAX = 20;

  CV = 0.0001;

  DV = 100000;

var

  N, A, B, Pa, C : array [1..4] of real;

  FU : array [1..8] of real;

  igo,i,j,ii:integer;

  k,base1,x,ipr,rp,t,ndelt,nsm0,kp,p:real;

  vm,v0,fo,df0,df1,df2,dx,rel,v,psum,calcrp,nsmr:real;

function func():real;

begin

  IF IGO > 2 THEN K := 1 ELSE K := -1;

  IF N[IGO] <= 0 then begin FUNC := 1; exit; end

  ELSE BASE1 := A[IGO] + K * N[IGO] * X;

  IF BASE1 > 0 THEN begin FUNC := deg(BASE1 , IPR); exit; end;

  func:=1;

end;

begin

  {//  ********************************************************

  //  *              Программа Р01 RAVN                      *

  //  ********************************************************

  //  *Расчет параметров состояния равновесия широкого класса*

  //  *            газофазных обратимых реакций              *

  //  *               aA + bB <=====> rR + sS                *

  //  ********************************************************

  // *          С ИСПОЛЬЗОВАНИЕМ МЕТОДА НЬЮТОНА             *

  // ********************************************************}

  writeln ('****************************************************');

  writeln ('*                РЕШЕНИЕ ЗАДАЧИ                    *');

  writeln ('****************************************************');  writeln;

  writeln ('************ВВОД ИСХОДНЫХ ДАННЫХ**************'); writeln;

  writeln ('ВВЕДИТЕ: '); writeln;

//210

  l1:

  write ('NA, NB, NR, NS - Стехиoметрические коэффициенты '); readln(N[1], N[2], N[3], N[4]);

  IF N[1] + N[2] = 0 THEN halt(1);

  writeln ('A0, B0, R0, S0 - исходные количества');

  write( '                 компонентов, Кмоль             '); readln(A[1], A[2], A[3], A[4]);

  write('KP - константа равновесия, МПа^[r+s-a-b]        '); readln(RP);

  write('T - температура, К                              '); readln(T);

  write( 'P - давление, МПа                               '); readln(P);

  write( 'X - нач. приближение характеристик равновесия   '); readln(X);

  { // Определение констант

   // ITMAX - Допустимое число операций

   // CV, DV - критерии сходимости}

  NDELT := N[3] + N[4] - (N[1] + N[2]);

  NSM0 := 0; FOR I := 1 TO 4 do NSM0 := NSM0 + A[I];

  VM := 22.4 * T * 0.1 / 273 / P; V0 := NSM0 * VM;

  KP := RP / deg((8.310001E-03 * T),NDELT) * deg(VM,NDELT);

  writeln; writeln( '****************РЕЗУЛЬТАТЫ РАСЧЕТА*************'); writeln;

  writeln('Последовательные приближения'); writeln;

  writeln('I              ','X              ', 'DX              ', 'F[X]  );

  writeln;

  {Итерационный цикл}

  FOR II := 1 TO ITMAX do begin

    FOR J := 1 TO 4 do begin

      IPR := N[J];

      IGO := J;

      FU[J] := FUNC;

    end;

    FOR J := 5 TO 8 do begin

      IPR := N[J - 4] - 1;

      IGO := J - 4;

      FU[J] := FUNC;

    end;

    NSMR := NSM0 + NDELT * X;

    FO := KP * deg(NSMR , NDELT) * FU[1] * FU[2] - FU[3] * FU[4];

    DF1 := -KP * deg(NSMR , NDELT) * (N[2] * N[2] * FU[1] * FU[6] + N[1] * N[1] * FU[2] * FU[5]);

    DF2 := KP * NDELT * NDELT * deg(NSMR , (NDELT - 1)) * FU[1] * FU[2];

    DF0 := DF1 + DF2 - N[4] * N[4] * FU[3] * FU[8] - N[3] * N[3] * FU[4] * FU[7];

    IF ABS(DF0)<=1E-25 then begin

      writeln( ' ФУНКЦИЯ DF[X] ПРИБЛИЖАЕТСЯ К 0');

      GOTO l1;

    end;

    DX := FO / DF0;

    write(II:20, X:20:5); writeln(DX:20:5, FO:20:5);

    REL := ABS(DX / X);

    IF REL <= CV then GOTO l2;

    IF REL >= DV then begin

      writeln( 'МЕТОД РАСХОДИТСЯ');

      GOTO l1;

    end;

    X := X - DX;

  end;

  writeln( ' НЕТ СХОДИМОСТИ ПРИ ITMAX:=', ITMAX);

  GOTO l1;

    {РАСЧЕТ РАВНОВЕСНЫХ КОЛИЧЕСТВ, РАВНОВЕСНЫХ ПАРЦИАЛЬНЫХ

     ДАВЛЕНИЙ КОМПОНЕНТОВ И КОНСТАНТЫ РАВНОВЕСИЯ}

  //810

  l2:

  FOR I := 1 TO 4 do IF I <= 2 THEN B[I] := A[I] - N[I] * X ELSE B[I] := A[I] + N[I] * X;

  V := VM * NSMR;

  FOR I := 1 TO 4 do begin

    Pa[I] := B[I] * 8.310001E-03 * T / V;

    C[I] := B[I] / NSMR;

  end;

  PSUM := 0;

  FOR I := 1 TO 4 do PSUM := PSUM + Pa[I];

  writeln; writeln('РАВНОВЕСНЫЕ КОЛИЧЕСТВА КОМПОНЕНТОВ'); writeln;

  writeln( '  AEQ:=', B[1]:8:5, '  BEQ:=', B[2]:8:5, '  REQ:=', B[3]:8:5, '  SEQ:=', B[4]:8:5, '[кмоль]');

  writeln ; writeln('ОБЪЕМ НАЧАЛЬНОЙ СМЕСИ V0:=', V0:10:5, '[м^3]');

  writeln('ОБЪЕМ РАВНОВЕСНОЙ СМЕСИ V:=', V:10:5, '[м^3]');

  writeln('РАВНОВЕСНЫЕ ПАРЦИАЛЬНЫЕ ДАВЛЕНИЯ КОМПОНЕНТОВ');

  writeln; writeln('PA:=', Pa[1]:10:5, 'PB:=', Pa[2]:10:5, '[МПа]');

  writeln('РАВНОВЕСНОЕ ОБЩЕЕ ДАВЛЕНИЕ P:=', PSUM:10:5, '[МПа]');

  writeln('РАВНОВЕСНЫЕ МОЛЯРНЫЕ ДОЛИ КОМПОНЕНТОВ');

  writeln('CA:=', C[1]:10:5, 'CB:=', C[2]:10:5, 'CR:=', C[3]:10:5, 'CS:=', C[4]:10:5);

  FOR I := 1 TO 4 do IF Pa[I] = 0 THEN Pa[I] := 1;

  CALCRP := deg(Pa[3] , N[3]) * deg(Pa[4] , N[4]) / (deg(Pa[1] , N[1]) * deg(Pa[2] , N[2]));

  writeln( 'РАСЧИТАННАЯ КР:=', CALCRP:10:5, '[МПа ^ [', NDELT, ']]');

  writeln ; writeln( 'ВЫПОЛНЕНИЕ ПРОГРАММЫ ЗАВЕРШЕНО');

end.

 

unit math;

interface

function IntPower(const Base: real; const Exponent: Integer): real;

function deg(const Base, Exponent: real): real;

implementation

function IntPower(const Base: real; const Exponent: Integer): real;

asm

        mov     ecx, eax

        cdq

        fld1                      { Result := 1 }

        xor     eax, edx

        sub     eax, edx          { eax := Abs(Exponent) }

        jz      @@3

        fld     Base

        jmp     @@2

@@1:    fmul    ST, ST            { X := Base * Base }

@@2:    shr     eax,1

        jnc     @@1

        fmul    ST(1),ST          { Result := Result * X }

        jnz     @@1

        fstp    st                { pop X from FPU stack }

        cmp     ecx, 0

        jge     @@3

        fld1

        fdivrp                    { Result := 1 / Result }

@@3:

        fwait

end;

function deg(const Base, Exponent: real): real;

begin

  if Exponent = 0.0 then

    Result := 1.0               { n**0 = 1 }

  else if (Base = 0.0) and (Exponent > 0.0) then

    Result := 0.0               { 0**n = 0, n > 0 }

  else if (Frac(Exponent) = 0.0) and (Abs(Exponent) <= MaxInt) then

    Result := IntPower(Base, Integer(Trunc(Exponent)))

  else

    Result := Exp(Exponent * Ln(Base))

end;

begin

end.

 

І.4. Програма на мові C++

#include <math.h>

#include <stdio.h>

#include <stdlib.h>

//********************************************************

//*              Программа Р01 RAVN                      *

//********************************************************

//*Расчет параметров состояния равновесия широкого класса*

//*            газофазных обратимых реакций              *

//*               aA + bB <=====> rR + sS                *

//********************************************************

//*          С ИСПОЛЬЗОВАНИЕМ МЕТОДА НЬЮТОНА             *

//********************************************************

double FUNC1050(double IPR, int IGO, double X, double *N, double *A)

{

      double K;

      if (IGO > 1) K = 1;

      else K = -1;

      double BASE1;

      if (N[IGO] > 0.0)

      {

             BASE1 = A[IGO] + K*N[IGO]*X;

             if (BASE1 > 0.0) return pow(BASE1,IPR);

      }

      return 1;

}

void main(void)

{

      double N[4],A[4],FU[8],B[4],P[4],C[4];

      double RP,T,_P,X;

      //печать сообщений об ошибке и возврат к началу программы

      printf("****************************************************\n");

      printf("*                РЕШЕНИЕ ЗАДАЧИ                    *\n");

      printf("****************************************************\n\n");

      printf("****************ВВОД ИСХОДНЫХ ДАННЫХ****************\n\n");

      printf("ВВЕДИТЕ:\n");

      while (true)

      {

      {

      float fN[4],fA[4];

      float fRP,fT,_fP,fX;

      printf("\nNA, NB, NR, NS - Стехиoметрические коэффициенты ");

      scanf("%f %f %f %f",&fN[0],&fN[1],&fN[2],&fN[3]);

      for (int i = 0; i < 4; i++) N[i] = fN[i];

      if (N[0] + N[1] == 0) return;

      printf("\nA0, B0, R0, S0 - исходные количества\n");

      printf("                 компонентов, Кмоль             ");

      scanf("%f %f %f %f",&fA[0],&fA[1],&fA[2],&fA[3]);

      for (i = 0; i < 4; i++) A[i] = fA[i];

      printf("\nKP - константа равновесия, МПа^(r+s-a-b)        "); scanf("%f",&fRP); RP = fRP;

      printf("\nT - температура, К                              "); scanf("%f",&fT); T = fT;

      printf("\nP - давление, МПа                               "); scanf("%f",&_fP); _P = _fP;

      printf("\nX - нач. приближение характеристик равновесия   "); scanf("%f",&fX); X = fX;

      }

      // Определение констант

      // ITMAX - Допустимое число операций

      // CV, DV - критерии сходимости

      int ITMAX = 20;

      double CV = 0.0001;

      double DV = 100000;

      double NDELT = (N[2] + N[3]) - (N[0] + N[1]);

      double NSM0 = 0;

      for (int i = 0; i < 4; i++) NSM0 = NSM0 + A[i];

      double VM = 22.4*T*0.1/273/_P;

      double V0 = NSM0*VM;

      double KP = RP/pow((8.310001E-03*T),NDELT)*pow(VM,NDELT);

      printf("\n");

      printf("****************РЕЗУЛЬТАТЫ РАСЧЕТА******************\n\n");

      printf("          Последовательные приближения\n\n");

      printf("  I         X               DX              F(X)\n\n");

      // Итерационный цикл

      for(int ii = 0; ii < ITMAX; ii++)

      {

             for (int j = 0; j < 4; j++)

                   FU[j] = FUNC1050(N[j],j,X,N,A);

             for (j = 4; j < 8; j++)

                   FU[j] =  FUNC1050(N[j-4]-1,j-4,X,N,A);

             double NSMR = NSM0 + NDELT * X;

             double FO = KP*pow(NSMR,NDELT)*FU[0]*FU[1]-FU[2]*FU[3];

             double DF1 = -KP*pow(NSMR,NDELT)*(N[1]*N[1]*FU[0]*FU[5]+N[0]*N[0]*FU[1]*FU[4]);

             double DF2 = KP*NDELT*NDELT*pow(NSMR,(NDELT-1))*FU[0]*FU[1];

             double DF0 = DF1+DF2-N[3]*N[3]*FU[2]*FU[7]-N[2]*N[2]*FU[3]*FU[6];

             if (fabs(DF0) <= 1E-25)

             {

                   printf(" ФУНКЦИЯ DF(X) ПРИБЛИЖАЕТСЯ К 0\n");

                   break;

             }

             double DX = FO/DF0;

             printf("%2i   %1.6E   %1.6E   %1.6E\n",ii,X,DX,FO);

             double REL = fabs(DX/X);

             if (REL <= CV)

             {

                   //расчет равновесных количеств, равновесных парциальных

                   //давлений компонентов и константы равновесия

                   for (i = 0; i < 4; i++)

                         if (i < 2) B[i] = A[i] - N[i]*X;

                         else B[i] = A[i] + N[i]*X;

                   double V = VM * NSMR;

                   for (i = 0; i < 4; i++)

                   {

                         P[i] = B[i] * 8.310001E-03 * T / V;

                         C[i] = B[i] / NSMR;

                   }

 

                   double PSUM = 0.0;

                   for (i = 0; i < 4; i++)

                         PSUM = PSUM + P[i];

                   printf("\n");

                   printf("          РАВНОВЕСНЫЕ КОЛИЧЕСТВА КОМПОНЕНТОВ\n\n");

                   printf("AEQ=%f  BEQ=%f REQ=%f SEQ=%f [кмоль]\n\n",B[0],B[1],B[2],B[3]);

                   printf("     ОБЪЕМ НАЧАЛЬНОЙ СМЕСИ V0=%f [м^3]\n\n",V0);

                   printf("     ОБЪЕМ РАВНОВЕСНОЙ СМЕСИ V= %f [м^3]\n\n",V);

                   printf("     РАВНОВЕСНЫЕ ПАРЦИАЛЬНЫЕ ДАВЛЕНИЯ КОМПОНЕНТОВ\n\n");

                   printf("PA=%f PB=%f [МПа]\n\n",P[0],P[1]);

                   printf("     РАВНОВЕСНОЕ ОБЩЕЕ ДАВЛЕНИЕ P=%f [МПа]\n\n",PSUM);

                   printf("     РАВНОВЕСНЫЕ МОЛЯРНЫЕ ДОЛИ КОМПОНЕНТОВ\n\n");

                   printf("CA=%f CB=%f CR=%f CS=%f\n\n",C[0],C[1],C[2],C[3]);

                   for (i = 0; i < 4; i++)

                         if (P[i] == 0.0) P[i] = 1.0;

                   double CALCRP = pow(P[2],N[2])*pow(P[3],N[3])/(pow(P[0],N[0])*pow(P[1],N[1]));

                   printf("РАСЧИТАННАЯ КР=%f[МПа^%.0f]\n\n",CALCRP,NDELT);

                   printf("ВЫПОЛНЕНИЕ ПРОГРАММЫ ЗАВЕРШЕНО\n\n");

                   exit(0);

             }

             if (REL >= DV)

             {

                   printf("МЕТОД РАСХОДИТСЯ\n\n");

                   break;

             }

             X = X - DX;

      }

      //Печать сообщений об ошибке и возврат к началу программы

      printf(" НЕТ СХОДИМОСТИ ПРИ ITMAX=%i\n",ITMAX);

      }

}

 

І.5. Програма на мові DELPHI

unit lab1;

interface

uses

  Windows, Messages, SysUtils, Classes, Graphics, Controls, Forms, Dialogs,

  StdCtrls, Menus, ExtCtrls;

type

  mas1 = array[1..4] of real;

  mas2 = array[1..4] of integer;

  TForm1 = class(TForm)

    Edit1: TEdit;

    Edit2: TEdit;

    Edit3: TEdit;

    Edit4: TEdit;

    Label1: TLabel;

    Label2: TLabel;

    Label3: TLabel;

    Label4: TLabel;

    Edit5: TEdit;

    Edit6: TEdit;

    Edit7: TEdit;

    Edit8: TEdit;

    Label5: TLabel;

    Label6: TLabel;

    Edit9: TEdit;

    Label7: TLabel;

    Edit10: TEdit;

    Label8: TLabel;

    Edit11: TEdit;

    Label9: TLabel;

    Label10: TLabel;

    Edit12: TEdit;

    Button1: TButton;

    Button2: TButton;

    Button3: TButton;

    Label11: TLabel;

    MainMenu1: TMainMenu;

    N1: TMenuItem;

    N2: TMenuItem;

    N3: TMenuItem;

    N4: TMenuItem;

    N5: TMenuItem;

    N6: TMenuItem;

    Edit13: TEdit;

    Edit14: TEdit;

    Edit15: TEdit;

    Edit16: TEdit;

    Label12: TLabel;

    Label13: TLabel;

    Label14: TLabel;

    Label15: TLabel;

    Label16: TLabel;

    Bevel1: TBevel;

    Bevel2: TBevel;

    Label17: TLabel;

    Label22: TLabel;

    Edit17: TEdit;

    Edit18: TEdit;

    Edit19: TEdit;

    Edit20: TEdit;

    Label18: TLabel;

    Label19: TLabel;

    Label23: TLabel;

    Label24: TLabel;

    Edit23: TEdit;

    Label25: TLabel;

    Edit24: TEdit;

    Edit25: TEdit;

    Edit26: TEdit;

    Edit27: TEdit;

    Label26: TLabel;

    Label27: TLabel;

    Label28: TLabel;

    Label29: TLabel;

    Label30: TLabel;

    Edit28: TEdit;

    Bevel3: TBevel;

    Bevel4: TBevel;

    Label20: TLabel;

    Label21: TLabel;

    Label31: TLabel;

    Label32: TLabel;

    Label33: TLabel;

    Label34: TLabel;

    Label35: TLabel;

    Label36: TLabel;

    Edit21: TEdit;

    Label37: TLabel;

    procedure Button3Click(Sender: TObject);

    procedure Button1Click(Sender: TObject);

    procedure EditKeyPress(Sender: TObject; var key: char);

    procedure Button2Click(Sender: TObject);

    procedure N6Click(Sender: TObject);

  private

    { Private declarations }

  public

    { Public declarations }

  end;

 

var

  Form1: TForm1;

implementation

{$R *.DFM}

function Power(Base,Exponent: Real): Real; {Vozv. v stepen}

var

 Tmp: Real;

begin

 Power:=0;

 if Int(Exponent/2)=Exponent/2 then Base:=Abs(Base);

 if Exponent=0 then begin

  Power:=1;

  Exit

 end;

 if Base=0 then begin

  power:=0;

  exit

 end;

 if Base>0 then begin

  Tmp:=Exp((Abs(Exponent))*Ln(Base));

  if Exponent>0 then begin

   power:=Tmp;

   exit

  end;

  if Exponent<0 then begin

   power:=1/Tmp;

   exit

  end;

 end;

 if Base<0 then begin

  Tmp :=Exp((Abs(Exponent))*Ln(Abs(Base)));

  if Exponent>0 then begin

   power:=-Tmp;

   exit

  end;

  if Exponent<0 then begin

   power:=-1/Tmp;

   exit

  end;

 end;

end;

function func1(igo,ipr:integer; X: real; A: mas1; N: mas2): real;

var

 base1: real;

 k: integer;

begin

 if igo>2 then k:=1 else k:=-1;

 if N[igo]>0 then

  begin

   base1:=A[igo]+k*N[igo]*x;

   if base1>0 then func1:=power(base1,ipr) else func1:=1;

  end

 else func1:=1;

end;

 

procedure Psiho(Edit1,Edit2,Edit3,Edit4,Edit5,Edit6,Edit7,Edit8,

 Edit9,Edit10,Edit11,Edit12,Edit13,Edit14,Edit15,Edit16,Edit17,

 Edit18,Edit19,Edit20,Edit21,Edit23,Edit24,Edit25,Edit26,Edit27,

 Edit28: TEdit; Label11: Tlabel);

const

 cv=0.0001;

 dv=1000000;

var

 i,j,ii,ipr,igo,itmax:integer;

 A,B,p,C: mas1;

 N: mas2;

 X,T,RP,Kp,Vo,V,dX,F0,DF0,DF1,DF2,Rel,calcrp,Ndelt,Nsm0,Nsmr,Vm,pp,psum:real;

 fu:array[1..8] of real;

 errormsg: string;

 done: boolean;

begin

 {Стехиометрические коэффициенты}

 N[1]:=StrToInt(Edit1.Text);

 N[2]:=StrToInt(Edit2.Text);

 N[3]:=StrToInt(Edit3.Text);

 N[4]:=StrToInt(Edit4.Text);

 itmax:=StrToInt(Edit21.Text);

 if (N[1]+N[2])<>0 then

 begin

  {Исходные количества компонентов, Кмоль}

  A[1]:=StrToFloat(Edit5.Text);

  A[2]:=StrToFloat(Edit6.Text);

  A[3]:=StrToFloat(Edit7.Text);

  A[4]:=StrToFloat(Edit8.Text);

  RP:=StrToFloat(Edit9.Text);{Кр - константа равновесия, МПа}

  T:=StrToFloat(Edit10.Text);{Т - температура, К}

  pp:=StrToFloat(Edit11.Text);{р - давление, МПа}

  X:=StrToFloat(Edit12.Text);{Х - нач. приближение}

  Ndelt:=N[3]+N[4]-(N[1]+N[2]);{разница стехиометрических коэффициентов}

  Nsm0:=0;

  For i:=1 to 4 do Nsm0:=Nsm0+A[i];

  Vm:=22.4*T*0.1/273/pp;

  Vo:=Nsm0*Vm;

  Kp:=RP/power(((8.310001E-03)*T),Ndelt)*power(Vm,Ndelt);

  ii:=0;

  Done:=false;

  Repeat

  ii:=ii+1;

   For j:=1 to 4 do

   begin

    ipr:=N[j];

    igo:=j;

    Fu[j]:=func1(igo,ipr,X,A,N);

   end;

   For j:=5 to 8 do

   begin

    ipr:=N[j-4]-1;

    igo:=j-4;

    Fu[j]:=func1(igo,ipr,X,A,N);

   end;

   Nsmr:=Nsm0+Ndelt*X;

   F0:=Kp*power(Nsmr,Ndelt)* Fu[1]*Fu[2]-Fu[3]*Fu[4];

   DF1:=-Kp*power(Nsmr,Ndelt)*(N[2]*N[2]*Fu[1]*Fu[6]+N[1]*N[1]*Fu[2]*Fu[5]);

   DF2:=Kp*Ndelt*Ndelt*power(Nsmr,(Ndelt-1))*Fu[1]*Fu[2];

   DF0:=DF1+DF2-N[4]*N[4]*Fu[3]*Fu[8]-N[3]*N[3]*Fu[4]*Fu[7];

   dX:=F0/DF0;

   Rel:=ABS(dX/X);

   X:=X-dX;

   if rel<=cv then

   begin

    done:=true;

    errormsg:='';

    {Расчёт равновесного состояния}

    For i:=1 to 4 do if i<=2 then b[i]:=A[i]-N[i]*X else B[i]:=A[i]+N[i]*X;

    V:=Vm*Nsmr;

    For i:=1 to 4 do

    begin

     p[i]:=B[i]*(8.310001E-03)*T/V;

     C[i]:=B[i]/Nsmr;

    end;

    psum:=0;

    For i:=1 to 4 do psum:=psum+p[i];

    {Вывод результатов}

    Edit13.Text:=FloatToStr(B[1]);

    Edit14.Text:=FloatToStr(B[2]);

    Edit15.Text:=FloatToStr(B[3]);

    Edit16.Text:=FloatToStr(B[4]);

    Edit17.Text:=FloatToStr(Vo);

    Edit18.Text:=FloatToStr(V);

    Edit19.Text:=FloatToStr(p[1]);

    Edit20.Text:=FloatToStr(p[2]);

    Edit23.Text:=FloatToStr(psum);

    Edit24.Text:=FloatToStr(C[1]);

    Edit25.Text:=FloatToStr(C[2]);

    Edit26.Text:=FloatToStr(C[3]);

    Edit27.Text:=FloatToStr(C[4]);

    For i:=1 to 4 do if p[i]=0 then p[i]:=1;

    calcrp:=power(p[3],N[3])*power(p[4],N[4])/power(p[1],N[1])/power(p[2],N[2]);

    Edit28.Text:=FloatToStr(calcrp);

   end;

   if rel>=dv then

   begin

    errormsg:='Метод расходится!';

    done:=true;

   end;

   if (ii=itmax) and (not done) then

    errormsg:='Нет сходимости при itmax='+IntToStr(itmax);

  until done or (ii=itmax);

 end

 else errormsg:='Задайте стехиометрические коэффициенты!';

 label11.Caption:=errormsg;

end;

procedure TForm1.Button3Click(Sender: TObject);

begin

 Form1.Close;

end;

procedure TForm1.Button1Click(Sender: TObject);

begin

 if (Edit1.Text='') or (Edit2.Text='') or (Edit3.Text='') or (Edit4.Text='')

 or (Edit5.Text='') or (Edit6.Text='') or (Edit7.Text='') or (Edit8.Text='')

 or (Edit9.Text='') or (Edit10.Text='') or (Edit11.Text='')

 or (Edit12.Text='') or (Edit21.Text='')

 then

  begin

   MessageDlg('Исходные данные введены не полностью!',mtError,[mbCancel],0);

  end

 else

  Psiho(Edit1,Edit2,Edit3,Edit4,Edit5,Edit6,Edit7,Edit8,

  Edit9,Edit10,Edit11,Edit12,Edit13,Edit14,Edit15,Edit16,Edit17,

  Edit18,Edit19,Edit20,Edit21,Edit23,Edit24,Edit25,Edit26,Edit27,

  Edit28,Label11);

end;

 

procedure TForm1.EditKeyPress(Sender: TObject; var key: char);

var

 buf: string[20];

 dcm: char;

 l: boolean;

begin

 dcm:=decimalseparator;

 buf:='';

 l:=false;

 if Sender=Edit1 then begin buf:=Edit1.Text; l:=true; end;

 if Sender=Edit2 then begin buf:=Edit2.Text; l:=true; end;

 if Sender=Edit3 then begin buf:=Edit3.Text; l:=true; end;

 if Sender=Edit4 then begin buf:=Edit4.Text; l:=true; end;

 if Sender=Edit5 then buf:=Edit5.Text;

 if Sender=Edit6 then buf:=Edit6.Text;

 if Sender=Edit7 then buf:=Edit7.Text;

 if Sender=Edit8 then buf:=Edit8.Text;

 if Sender=Edit9 then buf:=Edit9.Text;

 if Sender=Edit10 then buf:=Edit10.Text;

 if Sender=Edit11 then buf:=Edit11.Text;

 if Sender=Edit12 then buf:=Edit12.Text;

 if Sender=Edit21 then begin buf:=Edit21.Text; l:=true; end;

 case key of

  '0'..'9',chr(8):;

  '-': key:=chr(0);

  ',','.' :    begin

                if (pos(dcm,buf)<>0) or l then key:=chr(0)

                else

                if key<>dcm then key:=dcm;

               end;

  chr(13): begin

            if Sender=Edit1 then Edit2.SetFocus;

            if Sender=Edit2 then Edit3.SetFocus;

            if Sender=Edit3 then Edit4.SetFocus;

            if Sender=Edit4 then Edit5.SetFocus;

            if Sender=Edit5 then Edit6.SetFocus;

            if Sender=Edit6 then Edit7.SetFocus;

            if Sender=Edit7 then Edit8.SetFocus;

            if Sender=Edit8 then Edit9.SetFocus;

            if Sender=Edit9 then Edit10.SetFocus;

            if Sender=Edit10 then Edit11.SetFocus;

            if Sender=Edit11 then Edit12.SetFocus;

            if Sender=Edit12 then Edit21.SetFocus;

            if Sender=Edit21 then

            begin

             Button1.SetFocus;

             Button1Click(Sender);

            end;

           end;

  else key:=chr(0);

 end;

end;

 

procedure TForm1.Button2Click(Sender: TObject);

begin

 Edit1.Text:='';

 Edit2.Text:='';

 Edit3.Text:='';

 Edit4.Text:='';

 Edit5.Text:='';

 Edit6.Text:='';

 Edit7.Text:='';

 Edit8.Text:='';

 Edit9.Text:='';

 Edit10.Text:='';

 Edit11.Text:='';

 Edit12.Text:='';

 Edit13.Text:='';

 Edit14.Text:='';

 Edit15.Text:='';

 Edit16.Text:='';

 Edit17.Text:='';

 Edit18.Text:='';

 Edit19.Text:='';

 Edit20.Text:='';

 Edit23.Text:='';

 Edit24.Text:='';

 Edit25.Text:='';

 Edit26.Text:='';

 Edit27.Text:='';

 Edit21.Text:='';

end;

 

procedure TForm1.N6Click(Sender: TObject);

begin

 MessageDlg( 'Ravnogaz v1.03'+#13+#13

             + 'Программа расчёта равновесного состояния'+#13

             +'газофазных реакций'+#13+#13

             +'Автор:'+#13+#13

             +'  Студент Кафедры ХПСМ ИХФ НТУУ "КПИ"'+#13

             +'  Погодаев Константин Сергеевич'+#13+#13

             +'  группа ЛП-01' ,

             MtInformation,[mbOK],0);

end;

end.

 

Додаток ІІ. Програма розрахунку кінетики складних хімічних реакцій

IІ.1. Програма на мові BASIC

10 REM  *********************************************************************

20 REM  *               ПРОГРАММА P03 KINRK                                 *

30 REM  *********************************************************************

40 REM  *                                                                   *

50 REM  *    РЕШЕНИЕ СИСТЕМЫ КИНЕТИЧЕСКИХ УРАВНЕНИЙ ХИМИЧЕСКОЙ РЕАКЦИИ      *

60 REM  *                 СО СЛОЖНОЙ СХЕМОЙ ПРЕВРАЩЕНИЯ                     *

70 REM  *                                                                   *

80 REM  *                                 K1                                *

90 REM  *                        A + B --------> R + S                      *

100 REM *                              <-------                             *

110 REM *                                 K11                               *

120 REM *                                                                   *

130 REM *                                  K2                               *

140 REM *                        B + R --------> S + E                      *

150 REM *                                                                   *

160 REM *********************************************************************

170 REM *               МЕТОДОМ РУНГЕ-КУТТА ЧЕТВЕРТОГО ПОРЯДКА              *

180 REM *     ДИФФЕРЕНЦИАЛЬНЫЕ УРАВНЕНИЯ ОПИСАНЫ В ПОДПРОГРАММЕ FUNC        *

190 REM *********************************************************************

200 PRINT "*****************************************************************"

210 PRINT "*                                                               *"

220 PRINT "*                 Р Е Ш Е Н И Е  З А Д А Ч И                    *"

230 PRINT "*                                                               *"

240 PRINT "*****************************************************************"

250 PRINT : PRINT "********* ВВОД ИСХОДНЫХ ДАННЫХ ********"

260 PRINT : PRINT "     ВВЕДИТЕ :": PRINT

270 INPUT "N  -число уравнений                                          "; N

280 DIM Y(N), VEC(N), G(N)

290 INPUT "K1,K11,K2 - константы скоростей хим. реакций, л/(моль*с)"; K1, K11, K2

300 INPUT "X - НАЧАЛО ПРОВЕДЕНИЯ РЕАКЦИИ, c                             "; X

310 INPUT "XK - ПРОДОЛЖИТЕЛЬНОСТЬ РЕАКЦИИ, c                            "; XK

320 PRINT "НАЧАЛЬНАЯ КОНЦЕНТРАЦИЯ КОМПОНЕНТОВ РЕАКЦИИ, моль/л"

330 FOR J = 1 TO N: PRINT TAB(47); "Y0("; J; ")="; : INPUT Y(J): NEXT J

340 PRINT "    ЗАДАЙТЕ :  ": PRINT

350 INPUT "H - ШАГ ИНТЕГРИРОВАНИЯ                                       "; H

360 INPUT "Q - ШАГ ПЕЧАТИ РЕЗУЛЬТАТОВ                                   "; Q

380 PRINT : PRINT "********* РЕЗУЛЬТАТЫ РАСЧЕТА *********": PRINT

390 PRINT " ИЗМЕНЕНИЕ КОНЦЕНТРАЦИЙ КОМПОНЕНТОВ РЕАКЦИИ Y(J), моль/л "

400 PRINT "             ВО ВРЕМЕНИ   X,   c "

410 PRINT

420 GOSUB 540

430 M% = Q / H

440 GOSUB 580: REM обращение к п/п Ругне-Кутта

450 PRINT USING "###.##  "; X;

460 FOR J = 1 TO N: PRINT USING "####.####  "; Y(J); : NEXT J

470 PRINT

480 IF X <= XK THEN GOTO 440

500 PRINT

510 PRINT "ВЫПОЛНЕНИЕ ПРОГРАММЫ ЗАВЕРШЕНО"

520 END

540 PRINT "  X  "; : FOR J = 1 TO N: PRINT "   Y("; J; ")"; : NEXT J

550 PRINT

570 RETURN

580 REM ********************************************************************

590 REM *       ПОДПРОГРАММА РЕАЛИЗАЦИИ МЕТОДА РУНГЕ-КУТТА RUNKUT          *

600 REM ********************************************************************

610 DIM T1(N), T2(N), T3(N), T4(N)

620 FOR J = 1 TO M%

630 FOR I = 1 TO N

640 G(I) = Y(I)

650 NEXT I

660 GOSUB 800: REM обращение к п/п FUNC

670 FOR I = 1 TO N: T1(I) = H * VEC(I): Y(I) = G(I) + T1(I) / 2: NEXT I

680 GOSUB 800

690 FOR I = 1 TO N: T2(I) = H * VEC(I): Y(I) = G(I) + T2(I) / 2: NEXT I

700 GOSUB 800

710 FOR I = 1 TO N: T3(I) = H * VEC(I): Y(I) = G(I) + T3(I): NEXT I

720 GOSUB 800

730 FOR I = 1 TO N: T4(I) = H * VEC(I)

740 G(I) = G(I) + (T1(I) + 2 * (T2(I) + T3(I)) + T4(I)) / 6: NEXT I

750 FOR I = 1 TO N: Y(I) = G(I): NEXT I

760 X = X + H

770 NEXT J: ERASE T1, T2, T3, T4

780 RETURN

800 REM ********************************************************************

810 REM *    ПОДПРОГРАММА FUNC ВЫЧИСЛЕНИЯ ПРАВЫХ ЧАСТЕЙ ДИФ. УРАВНЕНИЙ     *

820 REM ********************************************************************

830 VEC(1) = -K1 * Y(1) * Y(2) + K11 * Y(3) * Y(4)

840 VEC(2) = -K1 * Y(1) * Y(2) + K11 * Y(3) * Y(4) - K2 * Y(2) * Y(3)

850 VEC(3) = K1 * Y(1) * Y(2) - K11 * Y(3) * Y(4) - K2 * Y(2) * Y(3)

860 VEC(4) = K1 * Y(1) * Y(2) - K11 * Y(3) * Y(4) + K2 * Y(2) * Y(3)

870 VEC(5) = K2 * Y(2) * Y(3)

880 RETURN

 

IІ.2. Програма на мові FORTRAN

      PROGRAM      matm2

      REAL K1,K11,K2

      DIMENSION Y(5), VEC(5), G(5), T1(5), T2(5), T3(5), T4(5)

      OPEN (1,file='REZUL.txt')

      WRITE(*,*) '****************************************************'

      WRITE(*,*) '*                RESHENIE ZADACHI                  *'

      WRITE(*,*) '****************************************************'

      WRITE(*,*) '************VVOD ICXODNUX DANNUX**************'

      WRITE(*,*)'VVEDITE:'

      WRITE(*,*)'N - chislo uravnenii'

      READ(*,*) N

      WRITE(*,*)'K1,K11,K2 - konstantu skorostei him. reakcii,L/(mol*s)'

      WRITE(*,*) 'K1'

      READ(*,*) K1

      WRITE(*,*) 'K11'

      READ(*,*) K11

      WRITE(*,*) 'K2'

      READ(*,*) K2

      WRITE(*,*) 'X -    nachalo provedenia reakcii,s'

      READ(*,*) X

      WRITE(*,*) 'XK - prodoljitelnost reakcii,s'

      READ(*,*) XK

      WRITE(*,*) 'nachalnaya koncentraciya komponentov reakcii, mol/l'

      DO J=1,N

      WRITE(*,*)'YO(',J,')='

      READ(*,*) Y(J)

      END DO

      WRITE(*,*) 'ZADAITE'

      WRITE(*,*) 'H -    shag integrirovaniya'

      READ(*,*) H

      WRITE(*,*) 'Q - shag pechati rezultatov'

      READ(*,*) Q

      WRITE(1,*) '****************РЕЗУЛЬТАТЫ РАСЧЁТА******************'

      WRITE(1,*)'изменение концентрации компонентов реакции Y(J),моль/л'

      WRITE(1,*) '                  во времени X, с'

!     WRITE (1,*) 'X='

!     DO J=1,N

      WRITE (1,*)  '      X         Y(1)        Y(2)       Y(3)      Y(4)

     *         Y(5)'                 

!     END DO

      M=Q/H

10    GOTO 30

40    WRITE(1,*)       

!     DO J=1,N

      WRITE (1,3) X,Y(1),Y(2),Y(3),Y(4),Y(5)

3     format (3x,F9.4,3x,f8.4,3x,f8.4,3x,f8.4,3x,f8.4,3x,f8.4)

!     END DO

      IF(X.LT.XK)  GOTO 10

      WRITE(1,*) 'выполнение программы завершено'

      STOP

30    CONTINUE

      DO J=0,M

      DO I=1,N

      G(I)=Y(I)

      END DO

      CALL FUN (VEC,N,K1,K2,K11,Y)

      DO I=1,N

      T1(I)=H*VEC(I)

      Y(I)=G(I)+T1(I)/2

      END DO

      CALL FUN (VEC,N,K1,K2,K11,Y)

      DO I=1,N

      T2(I)=H*VEC(I)

      Y(I)=G(I)+T2(I)/2

      END DO

      CALL FUN (VEC,N,K1,K2,K11,Y)

      DO I=1,N

      T3(I)=H*VEC(I)

      Y(I)=G(I)+T3(I)/2

      END DO

      CALL FUN (VEC,N,K1,K2,K11,Y)

      DO I=1,N

      T4(I)=H*VEC(I)

      G(I)=G(I)+(T1(I)+2*(T2(I)+T3(I))+T4(I))/6

      END DO

      DO I=1,N

      Y(I)=G(I)

      END DO

      X=X+H

      END DO

      GOTO 40

      END

 

c     Подпрограмма вычисления правых частей диф. уравнений

      SUBROUTINE FUN (VEC,N,K1,K2,K11,Y)

      DIMENSION VEC(N),Y(4)

      REAL K1,K2,K11

      VEC(1)=-K1*Y(1)*Y(2)+K11*Y(3)*Y(4)

      VEC(2)=-K1*Y(1)*Y(2)+K11*Y(3)*Y(4)-K2*Y(2)*Y(3)

      VEC(3)=K1*Y(1)*Y(2)-K11*Y(3)*Y(4)-K2*Y(2)*Y(3)

      VEC(4)=K1*Y(1)*Y(2)-K11*Y(3)*Y(4)+K2*Y(2)*Y(3)

      VEC(5)=K2*Y(2)*Y(3)

      END

 

IІ.3. Програма на мові PASCAL

const

  N=5; {'N  -число уравнений}

 

var

  Y, VEC, G : array [1..N] of real;

  k1,k11,k2:real;

  x,xk,h,q:real;

  j,m:integer;

 

procedure du();

begin

  {// ********************************************************************

    // *    ПОДПРОГРАММА FUNC ВЫЧИСЛЕНИЯ ПРАВЫХ ЧАСТЕЙ ДИФ. УРАВНЕНИЙ     *

    // ********************************************************************}

  VEC[1] := -K1 * Y[1] * Y[2] + K11 * Y[3] * Y[4];

  VEC[2] := -K1 * Y[1] * Y[2] + K11 * Y[3] * Y[4] - K2 * Y[2] * Y[3];

  VEC[3] := K1 * Y[1] * Y[2] - K11 * Y[3] * Y[4] - K2 * Y[2] * Y[3];

  VEC[4] := K1 * Y[1] * Y[2] - K11 * Y[3] * Y[4] + K2 * Y[2] * Y[3];

  VEC[5] := K2 * Y[2] * Y[3];

end;

 

procedure rk();

var

  T1, T2, T3, T4:array [1..n] of real;

  i,j:integer;

begin

  {rk -580 // ********************************************************************

  // *       ПОДПРОГРАММА РЕАЛИЗАЦИИ МЕТОДА РУНГЕ-КУТТА RUNKUT          *

  // ********************************************************************}

  FOR J := 1 TO M do begin

    FOR I := 1 TO N do G[I] := Y[I];

    du;// обращение к п/п FUNC

    FOR I := 1 TO N do begin

        T1[I] := H * VEC[I];

        Y[I] := G[I] + T1[I] / 2;

    end;

    du;

    FOR I := 1 TO N do begin

        T2[I] := H * VEC[I];

        Y[I] := G[I] + T2[I] / 2;

    end;

    du;

    fOR I := 1 TO N do begin

        T3[I] := H * VEC[I];

        Y[I] := G[I] + T3[I];

    end;

    du;

    FOR I := 1 TO N do begin

        T4[I] := H * VEC[I];

        G[I] := G[I] + (T1[I] + 2 * (T2[I] + T3[I]) + T4[I]) / 6;

    end;

    FOR I := 1 TO N do Y[I] := G[I];

    X := X + H;

  end;

end;

 

begin

  {//  *********************************************************************

  //  *               ПРОГРАММА P03 KINRK                                 *

  //  *********************************************************************

  writeln ('**************************************************************');

  writeln ('*                 Р Е Ш Е Н И Е  З А Д А Ч И                 *');

  writeln ('**************************************************************');

  writeln;

  writeln ('********* ВВОД ИСХОДНЫХ ДАННЫХ ********');

  writeln('     ВВЕДИТЕ :');

  write('K1,K11,K2 - константы скоростей хим. реакций, л/(моль*с) ');  readln(K1, K11, K2);

  write('X - НАЧАЛО ПРОВЕДЕНИЯ РЕАКЦИИ, c                             ');  readln(X);

  write('XK - ПРОДОЛЖИТЕЛЬНОСТЬ РЕАКЦИИ, c                            '); readln(XK);

  writeln ('НАЧАЛЬНАЯ КОНЦЕНТРАЦИЯ КОМПОНЕНТОВ РЕАКЦИИ, моль/л');

  FOR J := 1 TO N do begin

    writeln('Y0(', J, ')=');

    readln(Y[J]);

  end;

  writeln ('    ЗАДАЙТЕ :  '); writeln;

  write('H - ШАГ ИНТЕГРИРОВАНИЯ                                       '); readln(H);

  write('Q - ШАГ ПЕЧАТИ РЕЗУЛЬТАТОВ                                   '); readln(Q);

  writeln; writeln('********* РЕЗУЛЬТАТЫ РАСЧЕТА *********'); writeln;

  writeln(' ИЗМЕНЕНИЕ КОНЦЕНТРАЦИЙ КОМПОНЕНТОВ РЕАКЦИИ Y(J), моль/л ');

  writeln('             ВО ВРЕМЕНИ   X,   c ');  writeln;

  write('  X  ');

  FOR J := 1 TO N do write ('   Y(', J, ')');  writeln;

  M := round(Q / H);

  repeat

    rk ;  // обращение к п/п Ругне-Кутта

    write (X:5:2);

    FOR J := 1 TO N do write(Y[J]:8:4);  writeln;

  until X > XK;

  writeln;

  writeln('ВЫПОЛНЕНИЕ ПРОГРАММЫ ЗАВЕРШЕНО');

end.

 

IІ.4. Програма на мові C++

#include <math.h>

#include <iostream.h>

 

//    *********************************************************************

//    *               ПРОГРАММА P03 KINRK                                 *

//    *********************************************************************

void DefRightPartEq(double *Y, double K1, double K11, double K2, double *VEC)

{

      VEC[0] = -K1*Y[0]*Y[1] + K11*Y[2]*Y[3];

      VEC[1] = -K1*Y[0]*Y[1] + K11*Y[2]*Y[3] - K2*Y[1]*Y[2];

      VEC[2] =  K1*Y[0]*Y[1] - K11*Y[2]*Y[3] - K2*Y[1]*Y[2];

      VEC[3] =  K1*Y[0]*Y[1] - K11*Y[2]*Y[3] + K2*Y[1]*Y[2];

      VEC[4] =  K2*Y[1]*Y[2];

}

void main(void)

{

      cout <<      "*************************************************" << endl;

      cout <<      "*                РЕШЕНИЕ ЗАДАЧИ                 *" << endl;

      cout <<      "*************************************************" << endl;

      cout <<      "****************ВВОД ИСХОДНЫХ ДАННЫХ*************" << endl;

      cout <<      "ВВЕДИТЕ:" << endl;

      int N;

      cout <<      "N - число уравнений                                      ";

      cin >> N;

      double *Y;Y = new double[N];

      double K1,K11,K2;

      cout <<      "K1,K11,K2 - константы скоростей хим. реакций, л/(моль*с) ";

      cin >> K1 >> K11 >> K2;

      double X;

      cout <<      "X - НАЧАЛО ПРОВЕДЕНИЯ РЕАКЦИИ, c                         ";

      cin >> X;

      double XK;

      cout <<      "XK - ПРОДОЛЖИТЕЛЬНОСТЬ РЕАКЦИИ, c                        ";

      cin >> XK;

      cout << "НАЧАЛЬНАЯ КОНЦЕНТРАЦИЯ КОМПОНЕНТОВ РЕАКЦИИ, моль/л" << endl;

      int i;

      for (i = 0; i <    N; i++)

      {

             cout <<      "                                   Y0(" << i << ")=";

             cin >> Y[i];

      }

      cout <<      "ЗАДАЙТЕ :  " << endl;        

      double H;

      cout <<      "    H - ШАГ ИНТЕГРИРОВАНИЯ                               ";

      cin >> H;

      double Q;

      cout <<      "    Q - ШАГ ПЕЧАТИ РЕЗУЛЬТАТОВ                           ";

      cin >> Q;

      cout << endl;     

      cout <<      "********* РЕЗУЛЬТАТЫ РАСЧЕТА *********" << endl;

      cout << " ИЗМЕНЕНИЕ КОНЦЕНТРАЦИЙ КОМПОНЕНТОВ Y(J), моль/л" << endl;

      cout <<      "             ВО ВРЕМЕНИ   X,   c" << endl;

      cout <<      endl;

      cout << "  X    ";

      for (i = 0; i < N; i++)

             cout << "    Y(" << i << ")  ";

      cout << endl;           

      int M = Q/H;

      do {

      double *T1; T1 = new double[N];

      double *T2;  T2 = new double[N];

      double *T3;  T3 = new double[N];

      double *T4;  T4 = new double[N];

      double *VEC;VEC= new double[N]; for (int i = 0; i < N; i++) VEC[i] = 0;

      double *G;   G  = new double[N];

      for (int j = 0; j < M; j++)

      {

             for (i = 0; i < N; i++) G[i] = Y[i];

             DefRightPartEq(Y,K1,K11,K2,VEC);

             for (i = 0; i < N; i++) 

             {

                   T1[i] = H * VEC[i];

                   Y[i] = G[i] + T1[i]/2;

             }

             DefRightPartEq(Y,K1,K11,K2,VEC);

             for (i = 0; i < N; i++) 

             {

                   T2[i] = H * VEC[i];

                   Y[i] = G[i] + T2[i]/2;

             }

             DefRightPartEq(Y,K1,K11,K2,VEC);

             for (i = 0; i < N; i++) 

             {           

                   T3[i] = H * VEC[i];

                   Y[i] = G[i] + T3[i];

             }

             DefRightPartEq(Y,K1,K11,K2,VEC);

             for (i = 0; i < N; i++) 

             {           

                   T4[i] = H * VEC[i];

                   G[i] = G[i] + (T1[i] + 2*(T2[i] + T3[i]) + T4[i])/6;

             }

             for (i = 0; i < N; i++) Y[i] = G[i];

             X = X + H;

      }

      delete [] VEC;

      delete [] G;

      delete [] T4;

      delete [] T3;

      delete [] T2;

      delete [] T1;

      long oldFlags = cout.flags(ios::skipws | ios::fixed);

      int oldPrecision = cout.precision(2);

      cout << " " << X;

      cout.precision(4);

      for (i = 0; i < N; i++)

             cout << "    " << Y[i];

      cout << endl;

      cout.precision(oldPrecision);

      cout.flags(oldFlags);

      } while (X <= XK);

      cout << "ВЫПОЛНЕНИЕ ПРОГРАММЫ ЗАВЕРШЕНО";

      cout << endl;

      delete [] Y;

}

 

IІ.5. Програма на мові DELPHI

unit Unit1;

 

interface

uses

  Windows, Messages, SysUtils, Classes, Graphics, Controls, Forms, Dialogs,

  StdCtrls, ExtCtrls, Menus;

type

  mas1 = array[1..5] of real;

  mas2 = array[0..1000] of real;

  mas3 = array[0..1000,1..5] of real;

  TForm1 = class(TForm)

    Label1: TLabel;

    Edit1: TEdit;

    Edit2: TEdit;

    Edit3: TEdit;

    Label2: TLabel;

    Label3: TLabel;

    Label4: TLabel;

    Label5: TLabel;

    Edit4: TEdit;

    Label6: TLabel;

    Label7: TLabel;

    Edit5: TEdit;

    Label8: TLabel;

    Label9: TLabel;

    Edit6: TEdit;

    Edit7: TEdit;

    Edit8: TEdit;

    Edit9: TEdit;

    Edit10: TEdit;

    Label10: TLabel;

    Label11: TLabel;

    Label12: TLabel;

    Label13: TLabel;

    Label14: TLabel;

    Label15: TLabel;

    Label16: TLabel;

    Edit11: TEdit;

    Edit12: TEdit;

    Label17: TLabel;

    Label18: TLabel;

    MainMenu1: TMainMenu;

    N1: TMenuItem;

    N2: TMenuItem;

    N3: TMenuItem;

    N4: TMenuItem;

    N5: TMenuItem;

    Bevel1: TBevel;

    Memo1: TMemo;

    Button1: TButton;

    Button2: TButton;

    Button3: TButton;

    procedure Button1Click(Sender: TObject);

    procedure Button2Click(Sender: TObject);

    procedure Button3Click(Sender: TObject);

    procedure EditKeyPress(Sender: TObject; var key: char);

    procedure N5Click(Sender: TObject);

  private

    { Private declarations }

  public

    { Public declarations }

  end;

 

var

  Form1: TForm1;

implementation

{$R *.DFM}

procedure func(k1,k11,k2: real; var y,g,vec: mas1);

begin

 vec[1]:=-k1*y[1]*y[2]+k11*y[3]*y[4];

 vec[2]:=-k1*y[1]*y[2]+k11*y[3]*y[4]-k2*y[2]*y[3];

 vec[3]:=k1*y[1]*y[2]-k11*y[3]*y[4]-k2*y[2]*y[3];

 vec[4]:=k1*y[1]*y[2]-k11*y[3]*y[4]+k2*y[2]*y[3];

 vec[5]:=k2*y[2]*y[3];

end;

procedure roonkoot(k1,k11,k2,h:real; n,m:integer; Var x:real; Var y: mas1);

var

 vec,g,t1,t2,t3,t4: mas1;

 i,j: integer;

begin

 for j:=1 to m do

 begin

  for i:=1 to n do g[i]:=y[i];

  func(k1,k11,k2,y,g,vec);

  for i:=1 to n do

  begin

   t1[i]:=h*vec[i];

   y[i]:=g[i]+t1[i]/2;

  end;

  func(k1,k11,k2,y,g,vec);

  for i:=1 to n do

  begin

   t2[i]:=h*vec[i];

   y[i]:=g[i]+t2[i]/2;

  end;

  func(k1,k11,k2,y,g,vec);

  for i:=1 to n do

  begin

   t3[i]:=h*vec[i];

   y[i]:=g[i]+t3[i];

  end;

  func(k1,k11,k2,y,g,vec);

  for i:=1 to n do

  begin

   t4[i]:=h*vec[i];

   g[i]:=g[i]+(t1[i]+2*(t2[i]+t3[i])+t4[i])/6;

  end;

  x:=x+h;

 end;

end;

 

procedure pr_calc(k1,k11,k2,xn,xk,h,q: real; var y: mas1);

const

 n=5;

var

 m,iy: integer;

 x: real;

 s: string;

begin

 s:='x        Y1      Y2      Y3      Y4      Y5';

 Form1.Memo1.Lines.Add(s);

 x:=xn;

 s:=FloatToStrF(x,ffFixed,8,4)+'  ';

 for iy:=1 to 5 do s:=s+FloatToStrF(y[iy],ffFixed,8,4)+'  ';

 Form1.Memo1.Lines.Add(s);

 m:=round(q/h);

 repeat

  roonkoot(k1,k11,k2,h,n,m,x,y);

  s:=FloatToStrF(x,ffFixed,8,4)+'  ';

  for iy:=1 to 5 do s:=s+FloatToStrF(y[iy],ffFixed,8,4)+'  ';

  Form1.Memo1.Lines.Add(s);

 until x > (xk-h)

end;

 

procedure TForm1.EditKeyPress(Sender: TObject; var key: char);

var

 buf: string[20];

 dcm: char;

begin

 dcm:=decimalseparator;

 buf:='';

 if Sender=Edit1 then buf:=Edit1.Text;

 if Sender=Edit2 then buf:=Edit2.Text;

 if Sender=Edit3 then buf:=Edit3.Text;

 if Sender=Edit4 then buf:=Edit4.Text;

 if Sender=Edit5 then buf:=Edit5.Text;

 if Sender=Edit6 then buf:=Edit6.Text;

 if Sender=Edit7 then buf:=Edit7.Text;

 if Sender=Edit8 then buf:=Edit8.Text;

 if Sender=Edit9 then buf:=Edit9.Text;

 if Sender=Edit10 then buf:=Edit10.Text;

 if Sender=Edit11 then buf:=Edit11.Text;

 if Sender=Edit12 then buf:=Edit12.Text;

 case key of

  '0'..'9',chr(8):;

  '-': key:=chr(0);

  ',','.' :    begin

                if (pos(dcm,buf)<>0) then key:=chr(0)

                else

                if key<>dcm then key:=dcm;

               end;

  chr(13): begin

            if Sender=Edit1 then Edit2.SetFocus;

            if Sender=Edit2 then Edit3.SetFocus;

            if Sender=Edit3 then Edit4.SetFocus;

            if Sender=Edit4 then Edit5.SetFocus;

            if Sender=Edit5 then Edit6.SetFocus;

            if Sender=Edit6 then Edit7.SetFocus;

            if Sender=Edit7 then Edit8.SetFocus;

            if Sender=Edit8 then Edit9.SetFocus;

            if Sender=Edit9 then Edit10.SetFocus;

            if Sender=Edit10 then Edit11.SetFocus;

            if Sender=Edit11 then Edit12.SetFocus;

            if Sender=Edit12 then

            begin

             Button1.SetFocus;

             Button1Click(Sender);

            end;

           end;

  else key:=chr(0);

 end;

end;

 

procedure TForm1.Button2Click(Sender: TObject);

begin

 Edit1.Text:='';

 Edit2.Text:='';

 Edit3.Text:='';

 Edit4.Text:='';

 Edit5.Text:='';

 Edit6.Text:='';

 Edit7.Text:='';

 Edit8.Text:='';

 Edit9.Text:='';

 Edit10.Text:='';

 Edit11.Text:='';

 Edit12.Text:='';

 Memo1.Text:='';

end;

 

procedure TForm1.Button1Click(Sender: TObject);

var

 y: mas1;

 k1,k11,k2,xn,xk,h,q: real;

begin

 Memo1.Text:='';

 if (Edit1.Text='') or (Edit2.Text='') or (Edit3.Text='') or (Edit4.Text='')

 or (Edit5.Text='') or (Edit6.Text='') or (Edit7.Text='') or (Edit8.Text='')

 or (Edit9.Text='') or (Edit10.Text='') or (Edit11.Text='') or (Edit12.Text='')

 then

  begin

   MessageDlg('Исходные данные введены не полностью!',mtError,[mbCancel],0);

  end

 else

 begin

 k1:=StrToFloat(Edit1.Text);

 k11:=StrToFloat(Edit2.Text);

 k2:=StrToFloat(Edit3.Text);

 xn:=StrToFloat(Edit4.Text);

 xk:=StrToFloat(Edit5.Text);

 h:=StrToFloat(Edit11.Text);

 q:=StrToFloat(Edit12.Text);

 y[1]:=StrToFloat(Edit6.Text);

 y[2]:=StrToFloat(Edit7.Text);

 y[3]:=StrToFloat(Edit8.Text);

 y[4]:=StrToFloat(Edit9.Text);

 y[5]:=StrToFloat(Edit10.Text);

 pr_calc(k1,k11,k2,xn,xk,h,q,y);

 end;

end;

 

procedure TForm1.Button3Click(Sender: TObject);

begin

 Form1.Close;

end;

 

procedure TForm1.N5Click(Sender: TObject);

begin

MessageDlg( 'RoonKoot v1.01'+#13+#13

             + 'Программа расчёта изменения во времени концентраций'+#13

             +'реагентов с применинием метода Рунге-Кутта'+#13+#13

             +'Автор:'+#13+#13

             +'  Студент Кафедры ХПСМ ИХФ НТУУ "КПИ"'+#13

             +'  Погодаев Константин Сергеевич'+#13+#13

             +'  группа ЛП-01' ,

             MtInformation,[mbOK],0);

end;

 

end.

Додаток ІІI. Розрахунок каскада реакторів.

Програма на мові BASIC

 

30 REM  *********************************************************************

40 REM  *                     ПРОГРАММА P06  RISKAS                         *

50 REM  *********************************************************************

60 REM  *    РАСЧЕТ ЧИСЛА РЕАКТОРОВ ИДЕАЛЬНОГО СМЕШЕНИЯ В КАСКАДЕ,          *

70 REM  *   ОБЕСПЕЧИВАЮЩЕМ ДОСТИЖЕНИЕ ЗАДАННОЙ СТЕПЕНИ ПРЕВРАЩЕНИЯ ПРИ      *

80 REM  *     ПРОВЕДЕНИИ РЕАКЦИИ А ==> R ПРОИЗВОЛЬНОГО ПОРЯДКА N            *

90 REM  *********************************************************************

100 REM *       С ИСПОЛЬЗОВАНИЕМ ИТЕРАЦИОННОГО МЕТОДА НЬЮТОНА               *

110 REM *********************************************************************

120 REM

130 PRINT "*****************************************************************"

140 PRINT "*                                                               *"

150 PRINT "*                 Р Е Ш Е Н И Е  З А Д А Ч И                    *"

160 PRINT "*                                                               *"

170 PRINT "*****************************************************************"

180 PRINT

190 PRINT "********************* ВВОД ИСХОДНЫХ ДАННЫХ **********************"

200 PRINT

210 DIM CA(20), V(20), PV(20), XAK(20), R(J)

220 PRINT "ВВЕДИТЕ:": PRINT

230 INPUT "VR - ОБЪЕМ ОДНОГО РЕАКТОРА, м^3                              "; VR

240 INPUT "VL - ОБЪЕМНЫЙ РАСХОД ПОТОКА, м^3/мин                         "; VL

250 INPUT "N - ПОРЯДОК РЕАКЦИИ                                          "; N

260 INPUT "KT - КОНСТ. СКОРОСТИ РЕАКЦИИ, кмоль^(1-N)/(л^(1-N)*мин)      "; KT

270 INPUT "CA0 - КОНЦЕНТРАЦИЯ РЕАГЕНТА А ВО ВХОДНОМ ПОТОКЕ, кмоль/л     "; CA0

275 PRINT "ЗАДАЙТЕ:": PRINT

280 INPUT "XK - КОНЕЧНАЯ СТЕПЕНЬ ПРЕВРАЩЕНИЯ                            "; XK

290 INPUT "ITMAX - МАКСИМАЛЬНОЕ ЧИСЛО ИТЕРАЦИЙ                       "; ITMAX

291 INPUT "EPS - ТОЧНОСТЬ                                            "; EPS

300 PRINT

310 PRINT "*************** РЕЗУЛЬТАТЫ РАСЧЕТА ******************"

320 PRINT

330 PRINT "V - ДОЛЯ НЕПРОРЕГРИРОВАВШЕГО РЕАГЕНТА В РЕАКТОРЕ"

340 PRINT "CA - КОНЦЕНТРАЦИЯ РЕАГЕНТА НА ВЫХОДЕ ИЗ РЕАКТОРА, кмоль/л"

360 PRINT "XAK - СТЕПЕНЬ ПРЕВРАЩЕНИЯ РЕАГЕНТА В КАСКАДЕ": PRINT

400 J = 1: CA(0) = CA0: PV(0) = 1

405 TAU = VR / VL

410 V(J) = 1

420 R = KT * TAU * CA(J - 1) ^ (N - 1)

421 REM

422 REM ПОСЛЕДОВАТЕЛЬНЫЕ ИТЕРАЦИИ

423 REM

430 FOR I = 1 TO ITMAX

440 DF = N * R * V(J) ^ (N - 1) + 1

450 F = R * V(J) ^ N + V(J) - 1

460 S = F / DF

470 IF S <= EPS THEN GOTO 510

480 V(J) = V(J) - S

490 NEXT I

500 PRINT "НЕТ СХОДИМОСТИ ПРИ ITMAX="; ITMAX: GOTO 700

510 CA(J) = CA(J - 1) * V(J)

520 PV(J) = PV(J - 1) * V(J)

530 XAK(J) = 1 - PV(J)

540 REM

550 REM ПРОВЕРКА УСЛОВИЯ ОКОНЧАНИЯ РАСЧЕТА

560 REM

570 IF XAK(J) < XK THEN J = J + 1: GOTO 410

590 PRINT "------------------ V ----------- CA ------------ XAK ---"

600 PRINT

610 FOR I = 1 TO J

620 PRINT I; "РЕАКТОР"; : PRINT USING "   +#.####^^^^"; V(I); CA(I); XAK(I)

630 PRINT

640 NEXT I

650 PRINT "-----------------------------------------------------------------"

660 PRINT : PRINT SPC(10); : PRINT "НЕОБХОДИМОЕ ЧИСЛО РЕАТОРОВ В КАСКАДЕ = "; J

670 PRINT

680 PRINT "-----------------------------------------------------------------"

690 PRINT : PRINT "ВЫПОЛНЕНИЕ ПРОГРАММЫ ЗАВЕРШЕНО"

700 END