Оптимальная температура проведения реакции

function DATA
global xa0 A E R tau T_a T_b;
Обнуляем вектор А и вектор Е
A(2,1)=zeros; E(2,1)=zeros;
Концентрация реагента А в долях в реакции A-P-S (мольные доли)
xa0=1;
Время пребывания в реакторе (мин)
tau=10;
Параметры уравнения Аррениуса для первой реакции  
a) Предэскпоненциальный множитель (мин^(-1)):
A(1,1)=70;
б) Энергия активации (кал/моль):
E(1,1)=2500; 
Параметры уравнения Аррениуса для второй реакции  
 a) Предэскпоненциальный множитель (мин^(-1)):
A(2,1)=100;
б) Энергия активации (кал/моль):
E(2,1)=5000; 
Левая граница температурного интервала исследования (К);
T_a=350;
Правая граница температурного интервала исследования (К);
T_b=370;
Универсальная газовая постоянная (кал/моль/К):
R=1.9872;
end

Программный код файла model2_stat_T .m - расчет выходных концентраций продуктов из реактора
Значения вектора х функции model2_stat_T  зависят от параметра « T»
function x= model2_stat_T(T)
global xa0 A E R tau;
С помощью цикла «for» рассчитываем константы реакции k1, k2
for i=1:2
    k(i)=A(i,1)*exp(-E(i,1)/R/T);
end
Обнуляем матрицу a и вектор b
a=zeros(2,2);b=zeros(2,1);
Для решения системы уравнений используем метод обратной матрицы. Преобразуем систему уравнений: 
 
чтобы получить коэффициенты матрицы:
 
Задание коэффициентов матрицы А (коэффициенты перед  и  ) и 
вектора В (коэффициенты правой части системы уравнений, 
после знака «равно»)  .
a(1,1)=1+k(1)*tau;a(1,2)=-k(2)*tau;b(1)=xa0;
a(2,1)=tau*k(1);a(2,2)=(-1-tau*k(2)); b(2)=0;
Определение выходных параметров модели. А=х*В. х=А/В. х=А-1*В
x=inv(a)*b;
Оператор inv(а) - поиск обратной матрицы А
end

Программный код файла GLAV_model2_grafik.m - главная  управляющая программа
clc;
clear all;
close all;
global TT xa xp T_a T_b;
 i=0;
DATA;
Задаем значение коэффициента T (температура, К) в интервале от T_a до T_b с шагом 1
for T=T_a:1:T_b
     i=i+1;   
 Вызываем функцию x= model2_stat_T(T) 
x= model2_stat_T(T);
Массив значений коэффициента Т(К) в интервале от T_a до T_b с шагом 1
TT(i)=T;
Массив значений коэффициента ха (мол.дол) - концентрации продукта А на выходе из реактора в заданном интервале температур от T_a до T_b
xa(i)=x(1);
Массив значений коэффициента хр (мол.дол) - концентрации продукта Р на выходе из реактора в заданном интервале температур от T_a до T_b
xp(i)=x(2);
end
REPORT;

function REPORT;
global xa0 A E tau TT  xa xp T_a T_b R;
Выводим на экран необходимый текст.
disp('____________________________________________________________ ');
disp('ПРОГРАММА РАСЧЕТА ОПТИМАЛЬНОЙ ТЕМПЕРАТУРЫ В ИЗОТЕРМИЧЕСКОМ РЕАКТОРЕ ИДЕАЛЬНОГО ПЕРЕМЕШИВАНИЯ ');
disp('                                                                 РЕАКЦИЯ  :  A  =  P ' );
disp('____________________________________________________________ ');
disp('Программа включает следующие файлы: GLAV_model2_grafik.m+DATA.m+model2_stat_T+REPORT.m');
 disp('___________________________________________________________ ');
disp('ИСХОДНЫЕ ДАННЫЕ ');
disp('-------------------------------------------------');
Оператор num2str(xa0,'%10.2f') переводит числовые значения в строковые.
Для этого необходимо записать оператор num2str, далее в скобках указываем переменную, которую необходимо перевести в строковый формат «xa0», затем ставим запятую «,» либо пробел  и в одинарных кавычках через точку указываем число символов (%10) в строковой переменной и количество знаков после запятой(2f).
disp(['1.Концентрация реагента A на входе в реактор ( xa0 ) = ' num2str(xa0,'%10.2f') ' мольные доли']);
disp('--------------------------------------------------------');
disp(['2.Предэкспоненциальный множитель первой реакции( A(1)) = ' num2str(A(1),'%10.2f') ' мин^(-1)']);
disp(['3.Энергия активация первой реакции( E(1)) = ' num2str(E(1),'%10.2f') ' кал/моль']);
disp(['4.Предэкспоненциальный множитель второй реакции( A(2)) = ' num2str(A(2),'%10.2f') ' мин^(-1)']);
disp(['5.Энергия активация второй реакции( E(1)) = ' num2str(E(2),'%10.2f') ' кал/моль']);
disp(['6.Время пребывания в реакторе ( tau) = ' num2str(tau,'%10.2f') ' мин']);
disp(['7.Левая граница температурного интервала исследования (К)=' num2str(T_a,'%10.2f') ' K']);
disp(['8.Правая граница температурного интервала исследования (К)=' num2str(T_b,'%10.2f') ' K']);
disp(['9.Универсальная газовая постоянная (кал/моль/К)=' num2str(R,'%10.4f') ' кал/моль/К']);
disp('____________________________________________________________ ');
disp('---------------------------------------------');
disp ('РЕЗУЛЬТАТЫ РАСЧЕТОВ');
disp('------------------------------------------------');
disp('--------------------------------------------');
disp ('1.Концентрации продуктов на выходе из реактора');
disp('--------------------------------------------------');
disp('_____________________________')
disp(' T(K)             xa(мол.д.)       xp(мол.д.) ');  
disp('_____________________________')
i=0;
Выводим таблицу данных. Переменные T, xa, xp;
for i=1:length(TT)
        TTT=TT(i);xaa=xa(i);xpp=xp(i);
    disp(sprintf('%10.3f\t %10.5f\t %10.5f',TTT,xaa,xpp));
end
Построение графиков;
hfig(1)=figure;
hfig(2)=figure;
Рисунок 1. График зависимости выхода продукта A от температуры в реакторе.
figure(hfig(1));
Построение графика
plot(TT,xa);
Название графика
title('График зависимости выхода продукта A от температуры в реакторе');
Подписи осей графика.
xlabel('T - K');ylabel('Выходные концентрации : xa  - мольные доли');
Наложение сетки на график 
grid on;
Рисунок 2. График зависимости выхода продукта P от температуры в реакторе
figure(hfig(2));
Построение графика
plot(TT,xp);
Название графика
title('График зависимости выхода продукта P от температуры в реакторе');
Подписи осей графика.
xlabel('T - K');ylabel('Выходные концентрации : xp - мольные доли');
Наложение сетки на график 
grid on;
end