3.2 Выбор метода решения и программы

Математическое описание, данного ХТП представляет собой систему нелинейных алгебраических уравнений, аналитическое решение которой является достаточно сложной задачей. Для Решения такого математического описания применяется метод Ньютона-Рафсона. При этом выбранный метод последовательно применяют к каждому, из уравнений системы с контролем сходимости каждой переменной к корню с заданной погрешностью. Он реализуется следующим алгоритмом:

1. Для всех факторов xi математического описания Fi(x1, x2,…, xn)=0 задают начальные приближения xi0 (x10, x20,…,xn0).

2. На основе исходной системы уравнений формируют матрицу Якоби (Якобиан):

… )

…

………………… (

… )

При этом частные производные находят численно по формуле:

где Hi – малое приращение xi (Hi=ε(xi)),

ε – точность поиска решения (погрешность).

3. Решают следующую систему линейных алгебраических уравнений относительно приращений Δxi:

.=

Решение этой системы дает значения приращений Δxi (Δx1, Δx2, …, Δxn).

4. Вычисляют новые (уточненные) значения Δxi по формуле:

, (3.4)

где m – номер итерации.

5. Для всех Δxi должно выполняться условие ׀Δxi׀≤ε:

Для решения созданного математического описания используется программный пакет Matlab 2011.

IV. Решение задачи

Решение системы нелинейных алгебраических уравнений в среде программирования Matlab осуществляется с помощью процедуры fsolve, входящей в пакет Optimization Toolbox. Процедура fsolve решает уравнения вида

или системы n нелинейных алгебраических уравнений c n неизвестными

.

Для решения системы нелинейных алгебраических уравнений в Matlab создают 2 функции:

· Первая (SystNonlinear) содержит систему уравнений, в ней задают массив значений констант скоростей и начальных концентраций реагентов

· Вторая (SystNonlinear 1) – решает эту систему уравнений для заданного диапазона значений X (где X – время или др. переменная) и выводит решение в командном окне.

В окне WorkSpace создают m-файл, содержащий систему нелинейных алгебраических уравнений . М-файл можно создать из пункта меню File (File ® New ® Function). Концентрации веществ А, В, С и D – CA, CB, CC, CD обозначим как Y1, Y2, Y3, Y4).

1) Функция - SystNonlinear имеет вид:

function F=SystNonlinear(Y,X)

%Система нелинейных алгебраических уравнений

%Содержит пять уравнений с пятью неизвестными Y(1), Y(2), Y(3)

K=[4 0.5 5 0.7] %Массив констант скорости реакций

Y0=[7 0 0 0 0]%Массив исходных значений конц. реагентов

F(1)=Y0(1)-Y(1)+X*((-K(1)*Y(1)+K(2)*Y(2)^2+K(4)*Y(4)*Y(5));

F(2)=Y0(2)-Y(2)+X*(2*K(1)*Y(1)-2*K(2)*Y(2)^2-K(3)*Y(2)+K(4)*Y(4)*Y(5));

F(3)=Y0(3)-Y(3)+X*(2*K(3)*Y(2)-2*K(4)*Y(4)*Y(5));

F(4)=Y0(4)-Y(4)+X*(K(3)*Y(3)- K(4)*Y(4)*Y(5));

F(5)=Y0(5)-Y(5)+X*(2*K(1)*Y(1)-2*K(2)*Y(2)^2).

End

2) Функциия SystNonlinear1имеет вид:

function F=SystNonlinear1(Y,X)

%Решение системы нелинейных алгебр. Уравнений

%для диапазона значений X от 0 до 2 с шагом 0.1

X=0:0.1:2

M=length(X);

for i=1:M

c=X(i)

options=optimset('Display','on'); %Опция выходного отображения

[Y]=fsolve(@SystNonlinear,[1 0 0],options,c)%вызов оптимизатора

Z1(i)=Y(1)%Массив значений Y1 для диапазона значений X

Z2(i)=Y(2)%Массив значений Y2 для диапазона значений X

Z3(i)=Y(3)%Массив значений Y3 для диапазона значений X

Z4(i)=Y(4)%Массив значений Y4 для диапазона значений X

Z5(i)=Y(5)%Массив значений Y5 для диапазона значений X

end

В тексте функцииSystNonlinear1 для решения системы нелинейных алгебраических уравнений использованы следующие параметры вызова процедуры fsolve:

options=optimset('Display','on'); %Опция выходного отображения

[Y] = fsolve(@SystNonlinear,Y0,options,variable) % вызов оптимизатора.

[Y] – означает вывод значений Y при заданном значении X;

fsolve – вызов программы для численного решения системы нелинейных алгебраических уравнений; @SystNonlinear – ссылка на функцию, содержащую систему уравнений, options – используемые опции.

Для получения решения системы уравнений в командном окне набирают команду

>> SystNonlinear1,

Она запускает решение системы алгебраических уравнений, описанных в

тексте функцииSystNonlinear.

Решение системы выводится в конце командного окна в массивах Z1, Z2 Z3 Z4 иZ5, которые содержат массивы значений концентраций реагентов Y1, Y2 Y3 Y4и Y5 (при значениях X=0; 0,1; 0,2; 0,3;…; 2,0).

V. Результаты расчета.

5.1 Таблица результатов расчета