Теперь рассмотрим еще одну программу по той же задаче.

Для построения графика (рис.3.7.1)необходимо выбрать исходную величину-аргумент (будет отображаться по горизонтальной оси), диапазон и шаг ее изменения, а также величину-функцию, зависящую от аргумента (вертикальная ось). Построенный график покажет, как будет изменяться зависимая величина при изменении аргумента (при этом значения всех остальных исходных величин, которые были введены в окна ввода исходных величин будут оставаться неизменными) – рис. 3.7.2.

Рис.3.7.1. Ввод исходных данных и вычисление результатов.

Для построения трехмерного графика также надо выбрать необходимые данные. По умолчанию стоит ноль (0) в окошке ввода номера аргумента. Пока там ноль, программа будет строить двумерный график.

Рис.3.7.2. Построение графика (в данном случае открыт из файла).

Для построения зависимостей концентраций от времени, следует в окошке рядом с кнопкой Граф 1-5 ввести цифру от 1 до 4 и нажать кнопку Граф 1-5. Будет построен (рис.3.7.3) график зависимости концентрации какого-либо из реагентов от времени.

Если ввести цифру 5 и нажать кнопку Граф 1-5, то будет построен график их всех одновременно.

Рис.3.7.3. График зависимости концентрации реагентов от времени.

 

Функция VvodIshodnDannih описанная ниже, обеспечивает считывание значений исходных данных из соответствующих «окон» ввода в интерфейсе, показанном выше на рис.3.7.1. Тэги окон ввода edit1, …, edit11. Как видно, нумерация идет подряд.

%Программа расчета оптимального времени пребывания (tau - час) в%изотермическом реакторе идеального смешения

function VvodIshodnDannih(hObject, eventdata, handles)

global xa0 k1 k2 k3 tau_a tau_b ii aa bb pp NS xb0;

%Схемареакции

NS=str2double(get(handles.edit9,'String'));

%Концентрация реагента в долях в реакции (мольные доли)

xa0=str2double(get(handles.edit1,'String'));

%Константыскоростейреакций (час^(-1))

k1=str2double(get(handles.edit2,'String'));

k2=str2double(get(handles.edit3,'String'));

if NS==2

   k3=str2double(get(handles.edit10,'String'));

   xb0=str2double(get(handles.edit11,'String'));

else

   k3=0;

   xb0=0;

end%if

%Константысхемыреакции

aa=str2double(get(handles.edit4,'String'));

bb=str2double(get(handles.edit5,'String'));

pp=str2double(get(handles.edit6,'String'));

%Левая граница поиска оптимального времени пребывания в реакторе (час)

tau_a=str2double(get(handles.edit7,'String'));

%Правая граница поиска оптимального времени пребывания в реакторе (час)

tau_b=str2double(get(handles.edit8,'String'));

ii=0;

end

 

ФункцияVivodArgumentovNaEkran, описанная ниже, обеспечивает вывод на экран, переменные выводятся в окна с тэгамиedit1, …, edit11, то есть те же самые, которые служат для ввода исходной информации, и в соответствии с пояснениями (названиями вводимых параметров), написанными на панели интерфейса рядом с окнами ввода исходных значений. При этом в зависимости от номера схемы (1 или 2) выводятся надписи с кратким описанием схемы реакции.

function VivodArgumentovNaEkran(hObject, eventdata, handles)

global xa0 k1 k2 k3 tau_a tau_b ii aa bb pp NS xb0;

%vivod argumentov na ekran  

%Концентрация реагента А в долях в реакции A-P-S (мольные доли)

S=sprintf('%g',xa0);

set(handles.edit1,'String',S);

%Константы скоростей реакций (час^(-1))

S=sprintf('%g',k1);

set(handles.edit2,'String',S);

S=sprintf('%g',k2);

set(handles.edit3,'String',S);

%Константысхемыреакции

S=sprintf('%g',aa);

set(handles.edit4,'String',S);

S=sprintf('%g',bb);

set(handles.edit5,'String',S);

S=sprintf('%g',pp);

set(handles.edit6,'String',S);

%Левая граница поиска оптимального времени пребывания в реакторе (час)

S=sprintf('%g',tau_a);

set(handles.edit7,'String',S);

%Правая граница поиска оптимального времени пребывания в реакторе (час)

S=sprintf('%g',tau_b);

set(handles.edit8,'String',S);

ii=0;

if NS==1

   SSS1='aA -(k1)-> pP+bB -(k2)-> S';

   S=sprintf('%s',SSS1);

   set(handles.text126,'String',S);

end%if

if NS==2

   S=sprintf('%g',k3);

   set(handles.edit10,'String',S);

   S=sprintf('%g',xb0);

   set(handles.edit11,'String',S);

   SSS2='aA =(k1->,<-k2)= pP';

   SSS3='P+bB -(k3)-> sS';

   S=sprintf('%s',SSS2);

   set(handles.text127,'String',S);

   S=sprintf('%s',SSS3);

   set(handles.text128,'String',S);

end%if

%Схемареакции

S=sprintf('%g',NS);

set(handles.edit9,'String',S);

end

 

ФункцияVivodResultNaEkran выводит результаты вычислений в окна вывода на панели «Результаты вычислений» (рис.3.7.1) с тэгамиedit18,edit19,иedit22, …,edit25, то есть пронумерованными (за исключением 20-го и 21-го) подряд. Вывод происходит в окна в соответствии с подписями рядом с ними (рис.3.7.1) с названиями выводимых расчетных параметров (результатов вычислений). 

function VivodResultNaEkran(hObject, eventdata, handles)

global t_opt phi_p_max xa_last xb_last xp_last xs_last;

%vivod resultatov na ekran

%Оптимальные значения режимных параметров реактора

S=sprintf('%g',t_opt);

set(handles.edit18,'String',S);

S=sprintf('%g',phi_p_max);

set(handles.edit19,'String',S);

%Выходные концентрации после максимального времени

S=sprintf('%g',xa_last);

set(handles.edit22,'String',S);

S=sprintf('%g',xb_last);

set(handles.edit23,'String',S);

S=sprintf('%g',xp_last);

set(handles.edit24,'String',S);

S=sprintf('%g',xs_last);

set(handles.edit25,'String',S);

end

 

function x= model1_stat_tau(tau)

%Программный код файла model1_stat_tau.m - расчет выходных концентраций%продуктов из реактора %Программа расчета оптимального времени пребывания (tau - час) в %изотермическом реакторе идеального смешения со стехиометрической схемой %реакции

МЕТОД НЬЮТОНА-РАФСОНА РЕШЕНИЕ СИСТЕМЫ НЕЛИНЕЙНЫХ УРАВНЕНИЙ

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%схема №1

% aA -(k1)-> pP + bB -(k2)-> S

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

global xa0 k1 k2 aa bb pp xaNR xbNR xpNR ii;

 

eps=0.0001;

 

function f1=fun1(xVector, tau)

   xa=xVector(1);

   f1=xa0-xa-aa*tau*k1*(xa^aa);

end% fun1

 

function f2=fun2(xVector, tau)

   xa=xVector(1);

   xb=xVector(2);

   xp=xVector(3);

   f2=-xp+tau*pp*(k1*(xa^aa)-k2*(xp^pp)*(xb^bb));

end% fun2

 

function f3=fun3(xVector, tau)

   xa=xVector(1);

   xb=xVector(2);

   xp=xVector(3);

   f3=-xb+tau*bb*(k1*(xa^aa)-k2*(xp^pp)*(xb^bb));

end% fun3

 

function d11=dfun1_dxa(xVector, tau)

   xa=xVector(1);

   xb=xVector(2);

   xp=xVector(3);

   d11=-1-tau*k1*(aa^2)*(xa^(aa-1));

end% dfun1_dxa

 

function d21=dfun2_dxa(xVector, tau)

   xa=xVector(1);

   xb=xVector(2);

   xp=xVector(3);

   d21=(tau*pp*k1*aa*xa^(aa-1));

end% dfun2_dxa

 

function d22=dfun2_dxb(xVector, tau)

   xa=xVector(1);

   xb=xVector(2);

   xp=xVector(3);

   d22=-(tau*pp*k2*(xp^pp)*bb*(xb^(bb-1)));

end% dfun2_dxb

 

function d23=dfun2_dxp(xVector, tau)

   xa=xVector(1);

   xb=xVector(2);

   xp=xVector(3);

   d23=-1-(tau*(pp^2)*k2*(xb^bb)*(xp^(pp-1)));

end% dfun3_dxp

 

function d31=dfun3_dxa(xVector, tau)

   xa=xVector(1);

   xb=xVector(2);

   xp=xVector(3);

   d31=(tau*bb*k1*aa*(xa^(aa-1)));

end% dfun3_dxa

 

function d32=dfun3_dxb(xVector, tau)

   xa=xVector(1);

   xb=xVector(2);

   xp=xVector(3);

   d32=-1-(tau*(bb^2)*k2*(xp^pp)*(xb^(bb-1)));

end% dfun3_dxb

 

function d33=dfun3_dxp(xVector, tau)

   xa=xVector(1);

   xb=xVector(2);

   xp=xVector(3);

   d33=-(tau*bb*k2*(xb^bb)*pp*xp^(pp-1));

end% dfun3_dxp

 

%%% NEWTON - RAFSON %%%

 

if ii==0

   xVector=[xa0-0.01;0.01;0.01];

else

   xVector=[xaNR;xbNR;xpNR];

end% if

 

for i=1:100

 

   f1=fun1(xVector, tau);

   f2=fun2(xVector, tau);

   f3=fun3(xVector, tau);

   fVector=[f1;f2;f3];

 

   d11=dfun1_dxa(xVector, tau);

   d21=dfun2_dxa(xVector, tau);

   d22=dfun2_dxb(xVector, tau);

   d23=dfun2_dxp(xVector, tau);

   d31=dfun3_dxa(xVector, tau);

   d32=dfun3_dxb(xVector, tau);

   d33=dfun3_dxp(xVector, tau);

   Jac=[d11, 0, 0; d21, d22, d23; d31, d32, d33];

 

   JObr=Jac^(-1);

 

   hVector=JObr*fVector;

 

if norm(hVector)<eps

break;

end% if

 

 

   xVector=xVector-hVector;

 

end% for i

 

xaNR=xVector(1);

xbNR=xVector(2);

xpNR=xVector(3);

ii=ii+1;  

 

x=xVector;

 

end

 

function x= model2_stat_tau(tau)

%Программный код файла model2_stat_tau.m - расчет выходных концентраций %продуктов из реактора %Программа расчета оптимального времени пребывания (tau - час) в %изотермическом реакторе идеального смешения со стехиометрической схемой %реакции

МЕТОД НЬЮТОНА-РАФСОНА РЕШЕНИЕ СИСТЕМЫ НЕЛИНЕЙНЫХ УРАВНЕНИЙ

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%схема №2

% aA =(<-k2,k1->)= pP

% P + bB -(k3)-> sS

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

global xa0 xb0 k1 k2 k3 aa bb pp xaNR xbNR xpNR ii;

 

eps=0.0001;

 

function f1=fun1(xVector, tau)

   xa=xVector(1);

   xp=xVector(3);

   f1=xa0-xa+tau*aa*(k2*(xp^pp)-k1*(xa^(aa)));

end% fun1   

 

function f2=fun2(xVector, tau)

   xa=xVector(1);

   xb=xVector(2);

   xp=xVector(3);

f2=-xp+tau*(pp*k1*(xa^aa)-k2*pp*(xp^pp)-k3*xp*(xb^bb));

end% fun2

 

function f3=fun3(xVector, tau)

   xa=xVector(1);

   xb=xVector(2);

   xp=xVector(3);

   f3=xb0-xb-tau*bb*k3*xp*(xb^bb);

end% fun3

 

function d11=dfun1_dxa(xVector, tau)

   xa=xVector(1);

   d11=-1-tau*k1*(aa^2)*(xa^(aa-1));

end% dfun1_dxa

 

function d13=dfun1_dxp(xVector, tau)

   xa=xVector(1);

   xp=xVector(3);

   d13=tau*pp*k2*aa*(xp^(pp-1));

end% dfun1_dxp

 

function d21=dfun2_dxa(xVector, tau)

   xa=xVector(1);

   xb=xVector(2);

   xp=xVector(3);

   d21=tau*pp*k1*aa*(xa^(aa-1));

end% dfun2_dxa

 

function d22=dfun2_dxb(xVector, tau)

   xa=xVector(1);

   xb=xVector(2);

   xp=xVector(3);

   d22=-tau*bb*k3*xp*(xb^(bb-1));

end% dfun2_dxb

 

function d23=dfun2_dxp(xVector, tau)

   xa=xVector(1);

   xb=xVector(2);

   xp=xVector(3);

   d23=-1-tau*((pp^2)*k2*(xp^(pp-1))+k3*(xb^(bb)));

end% dfun2_dxp

 

function d32=dfun3_dxb(xVector, tau)

   xa=xVector(1);

   xb=xVector(2);

   xp=xVector(3);

   d32=-1-tau*(bb^2)*k3*xp*(xb^(bb-1));

end% dfun3_dxb

 

function d33=dfun3_dxp(xVector, tau)

   xa=xVector(1);

   xb=xVector(2);

   xp=xVector(3);

   d33=-tau*bb*k3*(xb^bb);

end% dfun3_dxp

 

 

%%% NEWTON - RAFSON %%%

 

if ii==0

   xVector=[xa0-0.01;xb0-0.01;0.01];

else

   xVector=[xaNR;xbNR;xpNR];

end% if

 

for i=1:100

 

   f1=fun1(xVector, tau);

   f2=fun2(xVector, tau);

   f3=fun3(xVector, tau);

   fVector=[f1;f2;f3];

 

   d11=dfun1_dxa(xVector, tau);

   d13=dfun1_dxp(xVector, tau);

   d21=dfun2_dxa(xVector, tau);

   d22=dfun2_dxb(xVector, tau);

   d23=dfun2_dxp(xVector, tau);

   d32=dfun3_dxb(xVector, tau);

   d33=dfun3_dxp(xVector, tau);

   Jac=[d11, 0, d13; d21, d22, d23; 0, d32, d33];

 

   JObr=Jac^(-1);

 

   hVector=JObr*fVector;

 

if norm(hVector)<eps

break;

end% if

 

 

   xVector=xVector-hVector;

end% for i

 

xaNR=xVector(1);

xbNR=xVector(2);

xpNR=xVector(3);

ii=ii+1;  

 

x=xVector;

 

end

 

Функцияpushbutton22_Callbackвызываетсякнопкой Граф 1-5 стэгомpushbutton22.Для построения зависимостей концентраций от времени, следует в окошке с тэгом edit26рядом с кнопкой Граф 1-5 ввести цифру от 1 до 4 и нажать кнопку Граф 1-5. Будет построен (рис.3.7.3) график зависимости концентрации какого-либо из реагентов от времени. Если ввести цифру 5 и нажать кнопку Граф 1-5, то будет построен график их всех одновременно.

После построения графика функцияpushbutton22_Callback выводит на экран (вызывает функции вывода, описанные выше) исходные данные и результаты расчетов.

function pushbutton22_Callback(hObject, eventdata, handles)

global t xa xb xp phi_p;

nGraf=str2double(get(handles.edit26,'String'));

if nGraf==1

   plot(t,xa,'k-');

title('концентрация А');

   xlabel('t (сек)');

   ylabel('концентрация мольн.доли.');

grid on;

end%if

if nGraf==2

   plot(t,xb,'r-');

title('концентрация B');

   xlabel('t (сек)');

   ylabel('концентрация мольн.доли.');

grid on;

end%if

if nGraf==3

   plot(t,xp,'b-');

title('концентрация P');

   xlabel('t (сек)');

   ylabel('концентрация мольн.доли.');

grid on;

end%if

if nGraf==4

   plot(t,phi_p,'g-');

title('параметр phi p выхода продкта P');

   xlabel('t (сек)');

   ylabel('концентрация мольн.доли.');

grid on;

end%if

if nGraf==5

for j=1:4

switch j

case 1

               ssr='-';

               xx=xa;

case 2

               ssr='--';

               xx=xb;

case 3

               ssr='-.';

               xx=xp;

otherwise

end%switch

if j==1

else

           hold on;

end%if

       plot(t,xx,ssr);

title('Выход продукта -A, --B и целевого продукта -.-.P');

       xlabel('t (сек)');

ylabel('Выход (мол.дол.) в изотермическом реакторе');

grid on; 

end% for j

end%if

VivodArgumentovNaEkran(hObject, eventdata, handles);

VivodResultNaEkran(hObject, eventdata, handles);

end

 

ФункцияReshitZadachu имеет центральное (наиглавнейшее) значение в ходе решения поставленной задачи (об отыскании результатов исходя из известных начальных данных). ФункцияReshitZadachu получает в качестве аргумента функции векторVectorArg который содержит все основные исходные данные, последовательно расположенные как элементы вектора. ФункцияReshitZadachu возвращает (ее возвращаемым значением является) вектор результатовVectorRes, который содержит вычисляемые результаты, которые в дальнейшем другими функциями будут выводиться в качестве результатов в окнах интерфейса (рис.3.7.1) и сохраняться в файл (ниже будет описано).

ФункцияReshitZadachu в зависимости от номера схемы (введенного в соответствующее окно) вызывает функцию решения по Схеме 1 (model1_stat_tau) или по Схеме 2 (model2_stat_tau). Также функцияReshitZadachu определяет оптимальные параметры.

function [VectorRes]=ReshitZadachu(VectorArg)

%Программа расчета оптимального времени пребывания (tau - час) в%изотермическом реакторе идеального смешения

global t xa xb xp phi_p;

global xa0 xb0 k1 k2 k3 aa bb pp tau_a tau_b xaNR xbNR xpNR ii;

xa0=VectorArg(1); k1=VectorArg(2); k2=VectorArg(3); aa=VectorArg(4); bb=VectorArg(5); pp=VectorArg(6); tau_a=VectorArg(7); tau_b=VectorArg(8); k3=VectorArg(9); xb0=VectorArg(10); NS=VectorArg(11);

ii=0;

disp('podgotovka dannih dlya grafikov koncentracij');

i=0;

for tau=tau_a:0.1:tau_b

i=i+1;  

if NS==1

   x= model1_stat_tau(tau);

end%if

if NS==2

   x= model2_stat_tau(tau);

end%if

t(i)=tau;

xa(i)=x(1);

xb(i)=x(2);

xp(i)=x(3);

phi_p(i)=xp(i)/xa0;

end%for tau

xa_last=xa(i);

xb_last=xb(i);

xp_last=xp(i);

[phi_p_max, k]=max(phi_p);

t_opt=t(k);

%%%%%%%%%%%%%%%%%%

VectorRes=[t_opt phi_p_max xa_last xb_last xp_last];

end

 

ФункцияReshitZadachuG аналогична описанной выше функции ReshitZadachu, однаков отличие от нее ReshitZadachuGпредназначена для многократного вызова при расчете точек, например, для трехмерного графика.

function [VectorRes]=ReshitZadachuG(VectorArg)

%Программа расчета оптимального времени пребывания (tau - час) в%изотермическом реакторе идеального смешения

global xa0 xb0 k1 k2 k3 aa bb pp tau_a tau_b xaNR xbNR xpNR ii;

ii=0;

xa0=VectorArg(1); k1=VectorArg(2); k2=VectorArg(3); aa=VectorArg(4); bb=VectorArg(5); pp=VectorArg(6); tau_a=VectorArg(7); tau_b=VectorArg(8); k3=VectorArg(9); xb0=VectorArg(10); NS=VectorArg(11);

disp('podgotovka dannih dlya grafikov koncentracij');

i=0;

for tau=tau_a:0.1:tau_b

i=i+1;  

if NS==1

   x= model1_stat_tau(tau);

end%if

if NS==2

   x= model2_stat_tau(tau);

end%if

t(i)=tau;

xa(i)=x(1);

xb(i)=x(2);

xp(i)=x(3);

phi_p(i)=xp(i)/xa0;

end%for tau

xa_last=xa(i);

xb_last=xb(i);

xp_last=xp(i);

[phi_p_max, k]=max(phi_p);

t_opt=t(k);

%%%%%%%%%%%%%%%%%%

VectorRes=[t_opt phi_p_max xa_last xb_last xp_last];

end

 

ФункцияPostrGraficстроитграфик. Каквидноизтестафункции, онаиспользуетимя (порядковыйномерподписанныйусоответствующегоокнавинтерфейсе) NameArgаргумента (отображаемогопогоризонтальнойоси), используетимя (порядковыйномерподписанныйусоответствующегоокнавинтерфейсе) NameVarзависимой переменной (отображаемойповертикальнойоси), шкалу по горизонтальной оси отзначенияLeftNameVar до значенияRightNameVar с шагомStepGr (шаг отображается в разметке горизонтальной оси и в шаге между точками графика).

Кроме того, функцияPostrGraficстроитграфик. Затем, функция подготавливает подписи к горизонтальной оси и к вертикальной оси и размещает их в соответствующих двух массивах. В каждом из массивов подписи из одинакового количества символов, включая пробелы. Последняя подпись в каждом массиве на один символ длиннее.

Затем, функция подготавливает подпись к горизонтальной оси, учитывая, что длина подписи принимается тринадцать символов (последняя подпись в массиве делается на один знак длиннее). Функция выбирает подпись - элемент массива с номером NameArg.

Также функция PostrGrafic подготавливает подпись к вертикальной оси, учитывая, что длина подписи принимается семнадцать символов (последняя подпись в массиве делается на один знак длиннее).Функция выбирает подпись - элемент массива с номером NameVar.

Затем функция PostrGrafic подписывает оси подготовленными для этого подписями.

КаквидноизтекстаArray_3D=0; zLlabel=0; NumberSteps_3D=0; LeftNameArg_3D=0; StepGr_3D=0; RightNameArg_3D=0;функция PostrGrafic присваивает значения равные нулю переменным, связанным с построением трехмерного графика. Это обеспечивает согласованную работу с некоторыми другими функциями (будут описаны ниже).

function PostrGrafic(hObject, eventdata, handles)

global NameVar LeftNameVar RightNameVar NameArg StepGr LeftNameArg_3D StepGr_3D RightNameArg_3D NameArg_3D;

global k1Array v1Array Array_3D xLlabel yLlabel zLlabel NumberSteps NumberSteps_3D;

global xa0 k1 k2 aa bb pp tau_a tau_b NS xb0 k3 t_opt phi_p_max xa_last xb_last xp_last;

VectorArg=[xa0 k1 k2 aa bb pp tau_a tau_b k3 xb0 NS];

VectorArg(NameArg)=LeftNameVar;

k1Array(1)=LeftNameVar;

[VectorRes]=ReshitZadachuG(VectorArg);

v1Array(1)=VectorRes(NameVar);

for i=2:NumberSteps+1

   k1Array(i)=k1Array(i-1)+StepGr;

   VectorArg(NameArg)=k1Array(i);       

   [VectorRes]=ReshitZadachuG(VectorArg);

   v1Array(i)=VectorRes(NameVar);

end

plot(k1Array,v1Array,'o-');

set(gca,'XGrid','on');

set(gca,'YGrid','on');

%podpisi k osyam

TextVectorArg=[' xa0(мол.д.) ',' k1(1/час) ',' k2(1/час) ',' aa     ',' bb     ',' pp     ',' t a(сек) ',' t b(сек) ',' k3(1/час) ',' xb0(мол.д.) '];

TextVectorRes=[' t opt(сек) ',' phi p max  ',' xa last(мол.д.) ',' xb last(мол.д.) ',' xp last(мол.д.) '];

xLlabel=' ';

for i=((NameArg-1)*13+1):(NameArg*13+1)

   xLlabel=strcat(xLlabel,TextVectorArg(i));

end

yLlabel=' ';

for i=((NameVar-1)*17+1):(NameVar*17+1)

   yLlabel=strcat(yLlabel,TextVectorRes(i));

end

xlabel(xLlabel);

ylabel(yLlabel);

Array_3D=0;

zLlabel=0;

NumberSteps_3D=0;

LeftNameArg_3D=0;

StepGr_3D=0;

RightNameArg_3D=0;

end

 

ФункцияPostrGrafic_3Dстроиттрехмерныйграфик. Функция использует для третьей оси (второй оси в горизонтальной плоскости) значение NameArg_3Dкак имя (номер указанный в подписях к элементам интерфейса) переменной по третьей оси (второй аргумент при построении графика зависимости от двух аргументов), шкалу по третьей оси от LeftNameArg_3D доRightNameArg_3D с шагом StepGr_3D. В конце работы функция присваиваетv1Array=0; для совместимости с некоторыми другими функциями (будут описаны ниже)

function PostrGrafic_3D(hObject, eventdata, handles)

global NameVar LeftNameVar RightNameVar NameArg StepGr LeftNameArg_3D StepGr_3D RightNameArg_3D NameArg_3D;

global k1Array v1Array Array_3D xLlabel yLlabel zLlabel NumberSteps NumberSteps_3D;

global xa0 k1 k2 aa bb pp tau_a tau_b NS xb0 k3 t_opt phi_p_max xa_last xb_last xp_last;

VectorArg=[xa0 k1 k2 aa bb pp tau_a tau_b k3 xb0 NS];

for i=1:NumberSteps+1

if i==1

       k1Array(i)=LeftNameVar; 

else

       k1Array(i)=k1Array(i-1)+StepGr;           

end%if

   VectorArg(NameArg)=k1Array(i);

for j=1:NumberSteps_3D+1

if j==1

           z1Array(j)=LeftNameArg_3D; 

else

           z1Array(j)=z1Array(j-1)+StepGr_3D;     

end%if

       VectorArg(NameArg_3D)=z1Array(j);

       [VectorRes]=ReshitZadachuG(VectorArg);

       Array_3D(j,i)=VectorRes(NameVar);

end%for j

end%for i

[x,y]=meshgrid(LeftNameVar:StepGr:RightNameVar,LeftNameArg_3D:StepGr_3D:RightNameArg_3D);

%x %числостолбцов=NumberSteps+1 числострок=NumberSteps_3D+1

%y %числостолбцов=NumberSteps+1 числострок=NumberSteps_3D+1

%Array_3D %числостолбцов=NumberSteps+1 числострок=NumberSteps_3D+1

meshc(x,y,Array_3D);

%podpisi k osyam

TextVectorArg=[' xa0(мол.д.) ',' k1(1/час) ',' k2(1/час) ',' aa     ',' bb     ',' pp     ',' t a(сек) ',' t b(сек) ',' k3(1/час) ',' xb0(мол.д.) '];

TextVectorRes=[' t opt(сек) ',' phi p max  ',' xa last(мол.д.) ',' xb last(мол.д.) ',' xp last(мол.д.) '];

xLlabel=' ';

for i=((NameArg-1)*13+1):(NameArg*13+1)

   xLlabel=strcat(xLlabel,TextVectorArg(i));

end

yLlabel=' ';

for i=((NameArg_3D-1)*13+1):(NameArg_3D*13+1)

   yLlabel=strcat(yLlabel,TextVectorArg(i));

end

zLlabel=' ';

for i=((NameVar-1)*17+1):(NameVar*17+1)

   zLlabel=strcat(zLlabel,TextVectorRes(i));

end

xlabel(xLlabel);

ylabel(yLlabel);

zlabel(zLlabel);

v1Array=0;

end

 

Функцияpushbutton13_Callback вызывается при нажатии кнопки интерфейса с надписью «ПОСТРОИТЬ ГРАФИК». Необходимо присвоить этой кнопке тэгpushbutton13 при создании интерфейса. Функция считывает значения переменных из окон интерфейса. NameVarиз окна «номер зависимой переменной» с тэгом edit31,LeftNameVarиз окна «минимальное значение аргумента» с тэгомedit29,RightNameVarиз окна «максимальное значение аргумента» с тэгомedit30,NameArgиз окна «номер аргумента» с тэгомedit28,StepGrиз окна«шаг аргумента» с тэгомedit32, а также с панели для построения трехмерного графика «Объемный график»считывает NameArg_3Dиз окна «номер аргумента» с тэгомedit33,LeftNameArg_3Dиз окна «минимальный» с тэгомedit34,RightNameArg_3Dиз окна «максимальный» с тэгомedit36,StepGr_3Dиз окна «шаг» с тэгомedit35. Необходимо присвоить этим окнам соответствующие тэги.

Если NameArg_3Dравеннулю, то вызывается функцияPostrGrafic для построения плоского графика. По умолчанию в окне панели для построения трехмерного графика «Объемный график» NameArg_3D«номер аргумента» (рис.3.7.1) стоит ноль. Его надо указать при создании интерфейса как значение ‘String’.

Если NameArg_3Dне равеннулю, то вызывается функцияPostrGrafic_3D для построения трехмерного графика.

function pushbutton13_Callback(hObject, eventdata, handles)

global PG;

if PG==1       

else

   clear global;

end%if

global NameVar LeftNameVar RightNameVar NameArg StepGr LeftNameArg_3D StepGr_3D RightNameArg_3D NameArg_3D;

global k1Array v1Array Array_3D xLlabel yLlabel zLlabel NumberSteps NumberSteps_3D;

global xa0 xb0 k1 k2 k3 aa bb pp tau_a tau_b NS t_opt phi_p_max xa_last xb_last xp_last;

NameVar=str2double(get(handles.edit31,'String'));

LeftNameVar=str2double(get(handles.edit29,'String'));

RightNameVar=str2double(get(handles.edit30,'String'));

NameArg=str2double(get(handles.edit28,'String'));

StepGr=str2double(get(handles.edit32,'String'));

VvodIshodnDannih(hObject, eventdata, handles);

NumberSteps=round((RightNameVar-LeftNameVar)/StepGr);

NameArg_3D=str2double(get(handles.edit33,'String'));

if NameArg_3D==0

   PostrGrafic(hObject, eventdata, handles);

else

LeftNameArg_3D=str2double(get(handles.edit34,'String'));

RightNameArg_3D=str2double(get(handles.edit36,'String'));

   StepGr_3D=str2double(get(handles.edit35,'String'));

   NumberSteps_3D=round((RightNameArg_3D-LeftNameArg_3D)/StepGr_3D);

   PostrGrafic_3D(hObject, eventdata, handles);

end%if

end

 

Функцияpushbutton3_Callbackвызываетсякнопкой «ВЫЧИСЛИТЬ». Необходимо этой кнопке присвоить тэгpushbutton3. Она обеспечивает вызов функцииVvodIshodnDannihто есть ввод исходных данных из окон интерфейса, вызов функцииReshitZadachu то есть вычисление, и вызов функцииVivodResultNaEkranто естьвывод результатов в окна интерфейса. В зависимости от номера схемы реакции, выводится соответствующая надпись на панель интерфейса.

function pushbutton3_Callback(hObject, eventdata, handles)

global t xa xb xp phi_p;

global t_opt phi_p_max xa_last xb_last xp_last;

global xa0 xb0 k1 k2 k3 aa bb pp tau_a tau_b NS;

VvodIshodnDannih(hObject, eventdata, handles);

VectorArg=[xa0 k1 k2 aa bb pp tau_a tau_b k3 xb0 NS];

[VectorRes]=ReshitZadachu(VectorArg);

t_opt=VectorRes(1); phi_p_max=VectorRes(2);

xa_last=VectorRes(3); xb_last=VectorRes(4);

xp_last=VectorRes(5);

VivodResultNaEkran(hObject, eventdata, handles);

if NS==1

   SSS1='aA -(k1)-> pP+bB -(k2)-> S';

   S=sprintf('%s',SSS1);

   set(handles.text126,'String',S);

end%if

if NS==2

   SSS2='aA =(k1->,<-k2)= pP';

   SSS3='P+bB -(k3)-> sS';

   S=sprintf('%s',SSS2);

   set(handles.text127,'String',S);

   S=sprintf('%s',SSS3);

   set(handles.text128,'String',S);

end%if

end

 

Функцияpushbutton2_Callback вызывается кнопкой «ВЫЧИСЛИТЬ И СОХРАНИТЬ», необходимо присвоить этой кнопке тэгpushbutton2, она обеспечивает вычисление результатов и сохранение в файл.

function pushbutton2_Callback(hObject, eventdata, handles)

global t xa xb xp phi_p;

global t_opt phi_p_max xa_last xb_last xp_last;

global xa0 xb0 k1 k2 k3 aa bb pp tau_a tau_b NS;

%Vichislit resultati

pushbutton3_Callback(hObject, eventdata, handles);

%otkrit fail dlya zapisi

[f,p]=uiputfile('C:\MATLAB701\work\OptimHTP_S1_7Res.txt','Окновыбора');

KudaZapisat=strcat(p,f);

f=fopen(KudaZapisat,'wt');

%gotovim dannije dlja avtomaticheskogo schitivanija

VectorArg=[xa0 k1 k2 aa bb pp tau_a tau_b k3 xb0 NS];

VectorRes=[t_opt phi_p_max xa_last xb_last xp_last];

for i=1:5

   fprintf(f,'%g\t',VectorRes(i));

end

fprintf(f,'\n');

for i=1:11

   fprintf(f,'%g\t',VectorArg(i));

end

fprintf(f,'\n');

lt=length(t);

fprintf(f,'%g\t',lt);

fprintf(f,'\n');

for i=1:lt

       fprintf(f,'%g\t',t(i));

end

fprintf(f,'\n');

for i=1:lt

       fprintf(f,'%g\t',xa(i));

end

fprintf(f,'\n');

for i=1:lt

       fprintf(f,'%g\t',xb(i));

end

fprintf(f,'\n');

for i=1:lt

       fprintf(f,'%g\t',xp(i));

end

fprintf(f,'\n');

for i=1:lt

       fprintf(f,'%g\t',phi_p(i));

end

fprintf(f,'\n');

%sohranit v fail resultati vichislenij

fprintf(f,'\nOptimHTP_S1_7 \n');

fprintf(f,'RESULTATI VICHISLENIY\n');

TextString=strcat(' t opt(сек) ',' phi p max  ',' xa last(мол.д.) ',' xb last(мол.д.) ',' xp last(мол.д.) ');

fprintf(f,TextString);

fprintf(f,'\n');

for i=1:5

   fprintf(f,'%g\t',VectorRes(i));

end

%sohranit v fail ishodnije argumenti

fprintf(f,'\nARGUMENTI VICHISLENIY\n');

TextString=strcat(' xa0(мол.д.) ',' k1(1/час) ',' k2(1/час) ',' aa     ',' bb     ',' pp     ',' t a(сек) ',' t b(сек) ',' k3(1/час) ',' xb0(мол.д.) ');

fprintf(f,TextString);

fprintf(f,'\n');

for i=1:11

   fprintf(f,'%g\t',VectorArg(i));

end

fclose(f);

end

 

Функция pushbutton14_Callbackвызывается кнопкой «ОТКРЫТЬ ИЗ ФАЙЛА», необходимо этой кнопке присвоить тэгpushbutton14, функция обеспечивает открытие из файла исходных данных и результатов и вывод их в окна интерфейса.

Графиконанестроит.

function pushbutton14_Callback(hObject, eventdata, handles)

clear global;

global t xa xb xp phi_p;

global t_opt phi_p_max xa_last xb_last xp_last;

global xa0 xb0 k1 k2 k3 aa bb pp tau_a tau_b NS;

%Otkrit is faila

[f,p]=uigetfile('C:\MATLAB701\work\OptimHTP_S1_7Res*.txt','Окновыбора');

OtkudaChitat=strcat(p,f);

f=fopen(OtkudaChitat,'rt');

%chitat is faila

for i=1:5

   VectorRes(i)=fscanf(f,'%g',1);

end

for i=1:11

   VectorArg(i)=fscanf(f,'%g',1);

end

lt=fscanf(f,'%g',1);

for i=1:lt

       t(i)=fscanf(f,'%g',1);

end

for i=1:lt

       xa(i)=fscanf(f,'%g',1);

end

for i=1:lt

       xb(i)=fscanf(f,'%g',1);

end

for i=1:lt

       xp(i)=fscanf(f,'%g',1);

end

for i=1:lt

       phi_p(i)=fscanf(f,'%g',1);

end

fclose(f);

t_opt=VectorRes(1); phi_p_max=VectorRes(2); xa_last=VectorRes(3); xb_last=VectorRes(4); xp_last=VectorRes(5);

xa0=VectorArg(1); k1=VectorArg(2); k2=VectorArg(3); aa=VectorArg(4); bb=VectorArg(5); pp=VectorArg(6); tau_a=VectorArg(7); tau_b=VectorArg(8); k3=VectorArg(9); xb0=VectorArg(10); NS=VectorArg(11);

%vivod na ekran

VivodArgumentovNaEkran(hObject, eventdata, handles);

VivodResultNaEkran(hObject, eventdata, handles);

end

 

Функцияpushbutton8_Callback вызывается при нажатии на кнопку «ПОСТРОИТЬ ГРАФИК И СОХРАНИТЬ», необходимо этой кнопке присвоить тэгpushbutton8, функция обеспечивает построение графика и его сохранение в файл.

function pushbutton8_Callback(hObject, eventdata, handles)

global NameVar LeftNameVar RightNameVar NameArg StepGr LeftNameArg_3D StepGr_3D RightNameArg_3D NameArg_3D;

global k1Array v1Array Array_3D xLlabel yLlabel zLlabel NumberSteps NumberSteps_3D;

global t_opt phi_p_max xa_last xb_last xp_last;

global xa0 xb0 k1 k2 k3 aa bb pp tau_a tau_b NS;

global PG;

%postroit grafic

PG=1;

pushbutton13_Callback(hObject, eventdata, handles);

PG=0;

%Sohranit grafic v fail

[f,p]=uiputfile('C:\MATLAB701\work\OptimHTP_S1_7Graf.txt','Окновыбора');

KudaZapisat=strcat(p,f);

f=fopen(KudaZapisat,'wt');

if Array_3D==0

%gotovim dannije dlja avtomaticheskogo chitivanija

   fprintf(f,'%d\t',NumberSteps);

for i=1:(NumberSteps+1)

       fprintf(f,'%g\t',k1Array(i));

end

   fprintf(f,'\n');

for i=1:(NumberSteps+1)

       fprintf(f,'%g\t',v1Array(i));

end

%podpisi k osyam grafica

   fprintf(f,'%s\t%s\n',xLlabel,yLlabel);

else

%gotovim dannije dlja avtomaticheskogo chitivanija

   NS=0;

   fprintf(f,'%d\n',NS);

fprintf(f,'%d\t%d\t%d\t%d\t%d\t%d\n',LeftNameVar,StepGr,RightNameVar,LeftNameArg_3D,StepGr_3D,RightNameArg_3D);

   fprintf(f,'\n');

for i=1:(NumberSteps_3D+1)

for j=1:(NumberSteps+1)           

          fprintf(f,'%g\t',Array_3D(i,j));

end%for j

       fprintf(f,'\n');

end

   fprintf(f,'\n');

%podpisi k osyam grafica

   fprintf(f,'%s\t%s\t%s\n',xLlabel,yLlabel,zLlabel);

end%if

%zapisivajem argumenti

VvodIshodnDannih(hObject, eventdata, handles);

VectorArg=[xa0 k1 k2 aa bb pp tau_a tau_b k3 xb0 NS];

for i=1:11

  fprintf(f,'%g\t',VectorArg(i));

end

fclose(f);

end

 

Функцияpushbutton15_Callback вызывается кнопкой «ОТКРЫТЬ ГРАФИК ИЗ ФАЙЛА», необходимо присвоить этой кнопке тэгpushbutton15, функция обеспечивает открытие графика из файла (то есть построение графика) и выводит в окна интерфейса значения исходных данных.

function pushbutton15_Callback(hObject, eventdata, handles)

clear global;

global t_opt phi_p_max xa_last xb_last xp_last;

global xa0 xb0 k1 k2 k3 aa bb pp tau_a tau_b NS;

%otkrit grafic is faila

[f,p]=uigetfile('C:\MATLAB701\work\OptimHTP_S1_7Graf*.txt','Окновыбора');

OtkudaChitat=strcat(p,f);

f=fopen(OtkudaChitat,'rt');

%chitat is faila

NumberSteps=fscanf(f,'%d',1);

if NumberSteps==0

   LeftNameVar=fscanf(f,'%g',1);

   StepGr=fscanf(f,'%g',1);

   RightNameVar=fscanf(f,'%g',1);

   LeftNameArg_3D=fscanf(f,'%g',1);

   StepGr_3D=fscanf(f,'%g',1);

   RightNameArg_3D=fscanf(f,'%g',1);

[x,y]=meshgrid(LeftNameVar:StepGr:RightNameVar,LeftNameArg_3D:StepGr_3D:RightNameArg_3D);

%числостолбцов=NumberSteps+1 числострок=NumberSteps_3D+1

   NumberSteps_3D=round((RightNameArg_3D-LeftNameArg_3D)/StepGr_3D);

   NumberSteps=round((RightNameVar-LeftNameVar)/StepGr);

for i=1:(NumberSteps_3D+1)

for j=1:(NumberSteps+1)           

          Array_3D(i,j)=fscanf(f,'%g',1);

end%for j

end% for i

%podpisi k osyam grafica

   xLlabel=fscanf(f,'%s',1);

   yLlabel=fscanf(f,'%s',1);

   zLlabel=fscanf(f,'%s',1);

%stroim grafic

   meshc(x,y,Array_3D);

   xlabel(xLlabel);

   ylabel(yLlabel);

   zlabel(zLlabel);

else

for i=1:NumberSteps+1

       k1Array(i)=fscanf(f,'%g',1);

end

for i=1:NumberSteps+1

       v1Array(i)=fscanf(f,'%g',1);

end

%podpisi k osyam grafica

   xLlabel=fscanf(f,'%s',1);

   yLlabel=fscanf(f,'%s',1);

%stroim grafic

   plot(k1Array,v1Array,'o-');

   set(gca,'XGrid','on');

   set(gca,'YGrid','on'); 

   xlabel(xLlabel);

   ylabel(yLlabel);

end%if

for i=1:11

  VectorArg(i)=fscanf(f,'%g',1);

end

fclose(f);

xa0=VectorArg(1); k1=VectorArg(2); k2=VectorArg(3);

aa=VectorArg(4); bb=VectorArg(5); pp=VectorArg(6); tau_a=VectorArg(7); tau_b=VectorArg(8);

k3=VectorArg(9); xb0=VectorArg(10); NS=VectorArg(11);

%vivod na ekrans);

VivodArgumentovNaEkran(hObject, eventdata, handles);

end

Для начинающих освоение МАТЛАБа, можно рекомендовать такую же по содержанию программу, но без визуального интерфейса.

Ниже изложено ее содержание.

Для наиболее начальной стадии изучения МАТЛАБа можно рекомендовать ее выполнить сначала не всю, а только следующие ее части:

ИСХОДНЫЕ ДАННЫЕ

ВЫЧИСЛИТЬ

ВЫВОД РЕЗУЛЬТАТОВ

Ниже будет пояснено, какие несколько функций для этого надо подготовить.

 

ИСХОДНЫЕ ДАННЫЕ

Функция DATA, необходимая для задания исходных данных, так как графического интерфейса с окнами ввода данных не предусмотрено. Данные помещаются в файлDATA, соответствующие переменные объявляются глобальными. В дальнейшем именно к файлуDATA будут обращаться различные функции, нуждающиеся в исходных данных, вместо прежнего, описанного выше, обращения к функции ввода данных из окон графического интерфейса.

function DATA

%Программный код файла DATA.m - задание исходной информации для расчетов;%Программа расчета оптимального времени пребывания (tau - час) в%изотермическом реакторе идеального смешения

global xa0 k1 k2 k3 tau_a tau_b ii aa bb pp NS xb0;

%Схема реакции  

NS=2;%1

%Концентрация реагента в долях в реакции (мольные доли)

xa0=0.213;

%Константы скоростей реакций (час^(-1))

k1=0.75;

k2=0.01;

   k3=0.333;

   xb0=0.92;

%Константы схемы реакции

aa=3;

bb=3;

pp=3;

%Левая граница поиска оптимального времени пребывания в реакторе (час)

tau_a=5;

%Правая граница поиска оптимального времени пребывания в реакторе (час)

tau_b=50;

ii=0;

end

 

ВЫЧИСЛИТЬ

ФункцияGLAV_Vichislit_pushbutton3_Callbackнесвязанани скакимиэлементами (кнопками) графического интерфейса, так как графического интерфейса нет, но выполняет те же самые задачи, что и раньше (как было описано в предыдущем варианте программы, где был графический интерфейс). Однако, вместо обращения к функции ввода данных из окон графического интерфейса, теперь функцияGLAV_Vichislit_pushbutton3_Callback обращается к функцииDATA, чтобы получить исходные данные. Объявленные в функцииGLAV_Vichislit_pushbutton3_Callback глобальные переменные в результате обращения к функцииDATA становятся содержащими те же значения, как одинаковые с ними по названию глобальные переменные в функцииDATA.

function GLAV_Vichislit_pushbutton3_Callback

global t xa xb xp phi_p;

global t_opt phi_p_max xa_last xb_last xp_last;

global xa0 xb0 k1 k2 k3 aa bb pp tau_a tau_b NS;

DATA;

VectorArg=[xa0 k1 k2 aa bb pp tau_a tau_b k3 xb0 NS];

[VectorRes]=ReshitZadachu(VectorArg);

t_opt=VectorRes(1); phi_p_max=VectorRes(2);

xa_last=VectorRes(3); xb_last=VectorRes(4);

xp_last=VectorRes(5);

VivodResultNaEkran;

end

 

function [VectorRes]=ReshitZadachu(VectorArg)

%Программа расчета оптимального времени пребывания (tau - час) в%изотермическом реакторе идеального смешения

global t xa xb xp phi_p;

global xa0 xb0 k1 k2 k3 aa bb pp tau_a tau_b xaNR xbNR xpNR ii;

xa0=VectorArg(1); k1=VectorArg(2); k2=VectorArg(3); aa=VectorArg(4);bb=VectorArg(5); pp=VectorArg(6); tau_a=VectorArg(7); tau_b=VectorArg(8); k3=VectorArg(9); xb0=VectorArg(10); NS=VectorArg(11);

ii=0;

disp('podgotovka dannih dlya grafikov koncentracij');

i=0;

for tau=tau_a:0.1:tau_b

i=i+1;  

if NS==1

   x= model1_stat_tau(tau);

end%if

if NS==2

   x= model2_stat_tau(tau);

end%if

t(i)=tau;

xa(i)=x(1);

xb(i)=x(2);

xp(i)=x(3);

phi_p(i)=xp(i)/xa0;

end%for tau

xa_last=xa(i);

xb_last=xb(i);

xp_last=xp(i);

[phi_p_max, k]=max(phi_p);

t_opt=t(k);

%%%%%%%%%%%%%%%%%%

VectorRes=[t_opt phi_p_max xa_last xb_last xp_last];

end

 

function x= model2_stat_tau(tau)

%Программный код файла model2_stat_tau.m - расчет выходных концентраций %продуктов из реактора %Программа расчета оптимального времени пребывания (tau - час) в %изотермическом реакторе идеального смешения со стехиометрической схемой %реакции

МЕТОД НЬЮТОНА-РАФСОНА СИСТЕМА НЕЛИНЕЙНЫХ УРАВНЕНИЙ

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%схема №2

% aA =(<-k2,k1->)= pP

% P + bB -(k3)-> sS

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

global xa0 xb0 k1 k2 k3 aa bb pp xaNR xbNR xpNR ii;

 

eps=0.0001;

 

function f1=fun1(xVector, tau)

   xa=xVector(1);

   xp=xVector(3);

   f1=xa0-xa+tau*aa*(k2*(xp^pp)-k1*(xa^(aa)));

end% fun1   

 

function f2=fun2(xVector, tau)

   xa=xVector(1);

   xb=xVector(2);

   xp=xVector(3);

f2=-xp+tau*(pp*k1*(xa^aa)-k2*pp*(xp^pp)-k3*xp*(xb^bb));

end% fun2

 

function f3=fun3(xVector, tau)

   xa=xVector(1);

   xb=xVector(2);

   xp=xVector(3);

   f3=xb0-xb-tau*bb*k3*xp*(xb^bb);

end% fun3

 

function d11=dfun1_dxa(xVector, tau)

   xa=xVector(1);

   d11=-1-tau*k1*(aa^2)*(xa^(aa-1));

end% dfun1_dxa

 

function d13=dfun1_dxp(xVector, tau)

   xa=xVector(1);

   xp=xVector(3);

   d13=tau*pp*k2*aa*(xp^(pp-1));

end% dfun1_dxp

 

function d21=dfun2_dxa(xVector, tau)

   xa=xVector(1);

   xb=xVector(2);

   xp=xVector(3);

   d21=tau*pp*k1*aa*(xa^(aa-1));

end% dfun2_dxa

 

function d22=dfun2_dxb(xVector, tau)

   xa=xVector(1);

   xb=xVector(2);

   xp=xVector(3);

   d22=-tau*bb*k3*xp*(xb^(bb-1));

end% dfun2_dxb

 

function d23=dfun2_dxp(xVector, tau)

   xa=xVector(1);

   xb=xVector(2);

   xp=xVector(3);

   d23=-1-tau*((pp^2)*k2*(xp^(pp-1))+k3*(xb^(bb)));

end% dfun2_dxp

 

function d32=dfun3_dxb(xVector, tau)

   xa=xVector(1);

   xb=xVector(2);

   xp=xVector(3);

   d32=-1-tau*(bb^2)*k3*xp*(xb^(bb-1));

end% dfun3_dxb

 

function d33=dfun3_dxp(xVector, tau)

   xa=xVector(1);

   xb=xVector(2);

   xp=xVector(3);

   d33=-tau*bb*k3*(xb^bb);

end% dfun3_dxp

 

%%% NEWTON - RAFSON %%%

 

if ii==0

   xVector=[xa0-0.01;xb0-0.01;0.01];

else

   xVector=[xaNR;xbNR;xpNR];

end% if

 

for i=1:100

 

   f1=fun1(xVector, tau);

   f2=fun2(xVector, tau);

   f3=fun3(xVector, tau);

   fVector=[f1;f2;f3];

 

   d11=dfun1_dxa(xVector, tau);

   d13=dfun1_dxp(xVector, tau);

   d21=dfun2_dxa(xVector, tau);

   d22=dfun2_dxb(xVector, tau);

   d23=dfun2_dxp(xVector, tau);

   d32=dfun3_dxb(xVector, tau);

   d33=dfun3_dxp(xVector, tau);

   Jac=[d11, 0, d13; d21, d22, d23; 0, d32, d33];

 

   JObr=Jac^(-1);

 

   hVector=JObr*fVector;

 

if norm(hVector)<eps

break;

end% if

 

   xVector=xVector-hVector;

end% for i

 

xaNR=xVector(1);

xbNR=xVector(2);

xpNR=xVector(3);

ii=ii+1;  

 

x=xVector;

 

end

 

function x= model1_stat_tau(tau)

%Программный код файла model1_stat_tau.m - расчет выходных концентраций%продуктов из реактора%Программа расчета оптимального времени пребывания (tau - час) в%изотермическом реакторе идеального смешения со стехиометрической схемой%реакции

МЕТОД НЬЮТОНА-РАФСОНА СИСТЕМА НЕЛИНЕЙНЫХ УРАВНЕНИЙ

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%схема №1

% aA -(k1)-> pP + bB -(k2)-> S

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

global xa0 k1 k2 aa bb pp xaNR xbNR xpNR ii;

 

eps=0.0001;

 

function f1=fun1(xVector, tau)

   xa=xVector(1);

   f1=xa0-xa-aa*tau*k1*(xa^aa);

end% fun1

 

function f2=fun2(xVector, tau)

   xa=xVector(1);

   xb=xVector(2);

   xp=xVector(3);

   f2=-xp+tau*pp*(k1*(xa^aa)-k2*(xp^pp)*(xb^bb));

end% fun2

 

function f3=fun3(xVector, tau)

   xa=xVector(1);

   xb=xVector(2);

   xp=xVector(3);

   f3=-xb+tau*bb*(k1*(xa^aa)-k2*(xp^pp)*(xb^bb));

end% fun3

 

function d11=dfun1_dxa(xVector, tau)

   xa=xVector(1);

   xb=xVector(2);

   xp=xVector(3);

   d11=-1-tau*k1*(aa^2)*(xa^(aa-1));

end% dfun1_dxa

 

function d21=dfun2_dxa(xVector, tau)

   xa=xVector(1);

   xb=xVector(2);

   xp=xVector(3);

   d21=(tau*pp*k1*aa*xa^(aa-1));

end% dfun2_dxa

 

function d22=dfun2_dxb(xVector, tau)

   xa=xVector(1);

   xb=xVector(2);

   xp=xVector(3);

   d22=-(tau*pp*k2*(xp^pp)*bb*(xb^(bb-1)));

end% dfun2_dxb

 

function d23=dfun2_dxp(xVector, tau)

   xa=xVector(1);

   xb=xVector(2);

   xp=xVector(3);

   d23=-1-(tau*(pp^2)*k2*(xb^bb)*(xp^(pp-1)));

end% dfun3_dxp

 

function d31=dfun3_dxa(xVector, tau)

   xa=xVector(1);

   xb=xVector(2);

   xp=xVector(3);

   d31=(tau*bb*k1*aa*(xa^(aa-1)));

end% dfun3_dxa

 

function d32=dfun3_dxb(xVector, tau)

   xa=xVector(1);

   xb=xVector(2);

   xp=xVector(3);

   d32=-1-(tau*(bb^2)*k2*(xp^pp)*(xb^(bb-1)));

end% dfun3_dxb

 

function d33=dfun3_dxp(xVector, tau)

   xa=xVector(1);

   xb=xVector(2);

   xp=xVector(3);

   d33=-(tau*bb*k2*(xb^bb)*pp*xp^(pp-1));

end% dfun3_dxp

 

%%% NEWTON - RAFSON %%%

 

%xVector=[xa;xb;xp];

if ii==0

   xVector=[xa0-0.01;0.01;0.01];

else

   xVector=[xaNR;xbNR;xpNR];

end% if

 

for i=1:100

 

   f1=fun1(xVector, tau);

   f2=fun2(xVector, tau);

   f3=fun3(xVector, tau);

   fVector=[f1;f2;f3];

 

   d11=dfun1_dxa(xVector, tau);

   d21=dfun2_dxa(xVector, tau);

   d22=dfun2_dxb(xVector, tau);

   d23=dfun2_dxp(xVector, tau);

   d31=dfun3_dxa(xVector, tau);

   d32=dfun3_dxb(xVector, tau);

   d33=dfun3_dxp(xVector, tau);

   Jac=[d11, 0, 0; d21, d22, d23; d31, d32, d33];

 

   JObr=Jac^(-1);

 

   hVector=JObr*fVector;

 

if norm(hVector)<eps

break;

end% if

 

   xVector=xVector-hVector;

 

end% for i

 

xaNR=xVector(1);

xbNR=xVector(2);

xpNR=xVector(3);

ii=ii+1;  

 

x=xVector;

 

end

 

ВЫВОДРЕЗУЛЬТАТОВ

ФункцияVivodResultNaEkranвыводитрезультаты вCommandWindow(так как графического интерфейса не предусмотрено).

function VivodResultNaEkran

global t xa xb xp phi_p NS;

%vivod resultatov na ekran

disp('ПРОГРАММА МОДЕЛИРОВАНИЯ ИЗОТЕРМИЧЕСКОГО РЕАКТОРА ИДЕАЛЬНОГО ПЕРЕМЕШИВАНИЯ ');

if NS==1

   disp('РЕАКЦИЯ: ');

   disp('aA -(k1)-> pP+bB -(k2)-> S');

end%if

if NS==2

   disp('РЕАКЦИЯ: ');

   disp('aA =(k1->,<-k2)= pP');

   disp('P+bB -(k3)-> sS');

end%if

disp ('РЕЗУЛЬТАТЫ РАСЧЕТОВ');

disp ('1.Концентрации продуктов на выходе из реактора');

disp(' tau(час) xa(мол.д.) xb(мол.д.) xp(мол.д.)'); 

i=0;

for i=1:length(t)

       tt=t(i);xaa=xa(i);xbb=xb(i); xpp=xp(i);

   disp(sprintf('%10.3f\t %10.5f\t %10.5f\t %10.5f',tt,xaa,xbb,xpp));

end

%Построение графиков;

hfig(1)=figure;

hfig(2)=figure;

hfig(3)=figure;

hfig(4)=figure;

figure(hfig(1));

plot(t,xa,'r-');

title('График зависимости выхода продукта A от времени пребывания в реакторе');

xlabel('tau - час');ylabel('Выходные концентрации: xa - мольные доли');

grid on;

figure(hfig(2));

plot(t,xb,'k-');