pfrtempopt_a.m

本代码为黄华江编著《实用化工计算机模拟—MATLAB在化学工程中的应用》的配套车程序

function PFRTempOpt_A
% 求反应管最优温度分布（PFR_ZhangJianhou4b.m)
%
%   Author: HUANG Huajiang
%   Copyright 2003 UNILAB Research Center, 
%   East China University of Science and Technology, Shanghai, PRC
%   $Revision: 1.0 $  $Date: 2003/02/14 $

global k k1 k2 k3

% 已知参数
k0 = 1; 
Z = 1;
epsilon = 0.02;     % 收敛相对精度
tspan1 = [0.00 0.14 0.26 0.38 0.52 0.64 0.76 0.88 1.00];
tspan2 = Inverse(tspan1);
k = [1.000 0.7 0.6 0.5 0.5 0.5 0.5 0.5 0.5];

m = 1;
kOld = 0;
while ( mean(abs((kOld-k)/k)) > epsilon )
    kOld = k;
    [t,x] = ode45(@CEquation,tspan1,[1.0 0.0],[],k,tspan1);    
    fprintf('\n')
    disp([tspan1' x])    
    k2 = Inverse(k);
     [t,rambda] = ode45(@rambdaEquation,tspan2,[0 1],[],k2,tspan2);
    fprintf('\n')
    disp([tspan2' rambda])
    rambda = Inverse(rambda)
    f = mean(abs(x(:,1).*(rambda(:,2)-rambda(:,1))+2*x(:,2).*k'.*(rambda(:,1)-3*rambda(:,2))))
    i = 2:length(k);
    k(i) = x(i,1).*(rambda(i,1)-rambda(i,2))./(2*x(i,2).*(rambda(i,1)-3*rambda(i,2)))
    m = m+1
end

% ------------------------------------------------------------------
function dxdt = CEquation(t,x,ki,tspan)                 % 浓度方程
k = spline(tspan,ki,t);
k1 = k;
k2 = k*k;
k3 = 2*k*k;
dxdt(1) = k2*x(2)-k1*x(1);
dxdt(2) = k1*x(1)-(k2+k3)*x(2);
dxdt = dxdt';

% ------------------------------------------------------------------
function f = rambdaEquation(t,rambda,k2,tspan)          % 伴随方程
k = spline(tspan,k2,t);
f(1) = k*(rambda(1)-rambda(2));
f(2) = k^2*(3*rambda(2)-rambda(1));
f = f';

% ------------------------------------------------------------------
function y = Inverse(x)                                 % 逆向排序
[row,col] = size(x);
if row==1
    i = col:-1:1;
    j = 1:col;
    y(j) = x(i);
elseif col==1
    i = row:-1:1;
    j = 1:row;
    y(j) = x(i);
else
    i = row:-1:1;
    j = 1:row;
    y(j,:) = x(i,:);
end