Решение систем нелинейных адгебраических уравнений методом Ньютона-Рафсона в среде MATLAB Решение нелинейных алгебраических уравнений и их систем в среде программирования MATLAB осуществляется с помощью процедуры fsolve, входящей в пакет Optimization Toolbox. Процедура fsolve решает уравнения вида  или системы n нелинейных алгебраических уравнений c n неизвестными  Чтобы получить профили концентраций реагентов в потоке ИС, нужно решить систему уравнений (24) при различных значениях времени t и получить для каждого момента времени концентрации веществ А, В, С и D – CA, CB, CC, CD (которые можно обозначить как YA, YB, YC, YD или Y1, Y2, Y3, Y4). Для решения системы нелинейных алгебраических уравнений в MATLAB создают 2 функции: первая (SystNonlinear) содержит систему уравнений, вторая (SystNonlinear1) – решает эту систему уравнений для заданного диапазона значений X (где X – время или др. переменная) и выводит решение в командном окне. В первой функции задают массив значений констант скоростей и начальных концентраций реагентов. Например, для реакции (15), протекающей в потоке ИС, текст функции SystNonlinear имеет вид: function F=SystNonlinear(Y,X) %Система нелинейных алгебраических уравнений %Содержит три уравнения с тремя неизвестными Y(1), Y(2), Y(3) K=[4 2 3 3 1.2] %Массив констант скорости реакций Y0=[5 2 0 0 0]%Массив исходных значений конц.реагентов F(1)=Y0(1)-Y(1)+X*(-K(1)*Y(2)+K(2)*Y(3)^0.7-K(3)*Y(1)^1.5); F(2)=Y0(2)-Y(2)+X*(-2*K(1)*Y(2)+2*K(2)*Y(3)^0.7); F(3)=Y0(3)-Y(3)+X*(3*K(1)*Y(2)-3*K(2)*Y(3)^0.7-K(4)*Y(3)*Y(4)+K(5)*Y(5)^2); F(4)=Y0(4)-Y(4)+X*(K(3)*Y(1)^1.5-K(4)*Y(3)*Y(4)+K(5)*Y(5)^2); F(5)=Y0(5)-Y(5)+X*(2*K(4)*Y(3)*Y(4)-2*K(5)*Y(5)^2); End Текст функции 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)%Массив значений переменной X присвоен с options=optimset('Display','on'); %Опция выходного отображения [Y]=fsolve(@SystNonlinear,[5 2 0 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)%Массив значений Y3 для диапазона значений X Z5(i)=Y(5)%Массив значений Y3 для диапазона значений X End End В тексте функцииSystNonlinear1 для решения системы нелинейных алгебраических уравнений использованы следующие параметры вызова процедуры fsolve: options=optimset('Display','on'); %Опция выходного отображения [Y] = fsolve(@SystNonlinear,Y0,options,variable) % вызов оптимизатора. [Y] –означает вывод значений Y при заданном значении X; fsolve – вызов программы для численного решения системы нелинейных алгебраических уравнений; @SystNonlinear –ссылка на функцию, содержащую систему уравнений, Y0 –массив исходных значений Y, необходимых для решения системы уравнений; options –используемые опции. variable –значения переменной X. Для получения решения системы уравнений в командном окне набирают команду >> SystNonlinear1, которая запускает решение системы алгебраических уравнений, описанных в тексте функцииSystNonlinear. Решение системы выводится в конце командного окна в массивах Z1, Z2, Z3, Z4иZ5которые содержат массивы значений концентраций реагентов Y1, Y2, Y3, Y4 и Y5 (при значениях X=0; 0,1; 0,2;…; 2,0), соответственно: Z1 = Columns 1 through 6 5.0000 3.0505 2.3555 1.9702 1.7176 1.5360 Columns 7 through 12 1.3976 1.2878 1.1980 1.1228 1.0588 1.0035 Columns 13 through 18 0.9551 0.9123 0.8741 0.8398 0.8087 0.7804 Columns 19 through 21 0.7545 0.7308 0.7088 Z2 = Columns 1 through 6 2.0000 1.2978 1.0491 0.9182 0.8375 0.7828 Columns 7 through 12 0.7432 0.7132 0.6897 0.6707 0.6550 0.6419 Columns 13 through 18 0.6307 0.6211 0.6127 0.6053 0.5988 0.5930 Columns 19 through 21 0.5877 0.5830 0.5787 Z3 = Columns 1 through 6 0 0.7795 0.8733 0.9081 0.9272 0.9391 Columns 7 through 12 0.9472 0.9530 0.9572 0.9603 0.9628 0.9646 Columns 13 through 18 0.9661 0.9673 0.9682 0.9690 0.9696 0.9701 Columns 19 through 21 0.9705 0.9708 0.9711 Z4 = Columns 1 through 6 0 1.3246 1.6160 1.7742 1.8846 1.9688 Columns 7 through 12 2.0361 2.0916 2.1385 2.1788 2.2139 2.2449 Columns 13 through 18 2.2725 2.2972 2.3196 2.3399 2.3585 2.3756 Columns 19 through 21 2.3914 2.4060 2.4196 Z5 = Columns 1 through 6 0 0.5476 1.1062 1.4293 1.6332 1.7733 Columns 7 through 12 1.8759 1.9544 2.0167 2.0673 2.1094 2.1451 Columns 13 through 18 2.1756 2.2021 2.2254 2.2460 2.2644 2.2809 Columns 19 through 21 2.2958 2.3094 2.3218 Полученные решения системы уравнений переносят в таблицу, приведенную выше, и строят графики зависимости концентраций реагентов от времени. |