fsolve regression

[John W. Eaton-3]

With Octave 3.0.5, fsolve computes a good solution for the following
problem.  The norm of the function values is approximately 1e-5.  With
3.2.3 and Octave from the current hg archive, the norm of the function
values is about 3500.

jwe



global  npts A Aint Bint Rint caf cbf ccf ...
        k10 E1 k20 E2 Ka0 Ea Kc0 Ec Da Db Dc T ...
        Ra Rb Rc Da Db Dc kma kmb kmc dcadr dcbdr dccdr ...
       
Rg  = 8.314;  % J/K mol
Mair= 28.96; % g/mol, mol weight air
%Pin   = 1.5*1.013e5;  % N/m^2
Pin   = 2.0*1.013e5;  % N/m^2
Pext  = 1.013e5; % minimum allowed exit pressure
Tin   = 550; % K
Rt  = 10/2; %cm, tube radius
At  = pi*Rt*Rt; % cm^2, tube cross-sectional area
Ta  = 325; % K, ambient temperature
U    = 5.5e-3; % cal/(cm^2 K s), heat transfer coefficient
delH1 = (-67.63e3); % cal/mol CO;
delH2 = (-460.4e3); % cal/mol C_3H_6;
Cp  = 0.25; % cal/(g K), heat capacity of air
vis = 0.028e-2; % g/(cm s), viscosity of air
bt  = 1; % bed to tube area ratio
%u   = 500/bt;  % cm/sec, feed gas velocity entering bed
u  = 75/bt;  % cm/sec, feed gas velocity entering bed
Ac  = bt*At;  % cm^2 area of bed, 4 x area of tube
Qin = u*Ac;  % cm^3/sec, feed volumetric flowrate
P   = Pin;
T   = Tin;
Q   = Qin;
cfin= P/(Rg*T)*1e-6; % mol/cm^3
Rp  = 0.175; % cm radius of catalyst particle
a   = Rp/3;
rhob = 0.51; % g/cm^3 bed density
rhop = 0.68; % g/cm^3 particle density
epsb = 1 - rhob/rhop; % bed porosity, dimensionless
%epsb = 0.4;
xc  = 0.996;

cafin = cfin*0.02;
cbfin = cfin*0.03;
ccfin = cfin*5e-4;

Nafin = Q*cafin;
Nbfin = Q*cbfin;
Ncfin = Q*ccfin;
Nfin  = [Nafin; Nbfin; Ncfin];
Ntfin = Q*cfin;
massf = Ntfin*Mair; % mass flowrate, g/s, remains constant
Cptot = massf*Cp;

alpha = 1;
k10  = alpha*6.802e16*2.6e6*80/100*0.05/100;  %mol/cm^3 s
k20  = alpha*1.416e18*2.6e6*80/100*0.05/100;  %mol/cm^3 s
E1   = 13108; %K
E2   = 15109; %K
Ka0  = 8.099e6; % cm^3/mol
Kc0  = 2.579e8; % cm^3/mol
Ea   = - 409; %K
Ec   = 191; %K
Da   = 0.0487; % cm^2/s
Db   = 0.0469; % cm^2/s
Dc   = 0.0487; % cm^2/s
kma  = 0.4*9.76;   % cm/s
kmb  = 0.4*10.18;  % cm/s
kmc  = 0.4*9.76;   % cm/s

%%
%% calculate initial reaction rates
%%
k1      = k10*exp(-E1/Tin);
k2      = k20*exp(-E2/Tin);
Ka      = Ka0*exp(-Ea/Tin);
Kc      = Kc0*exp(-Ec/Tin);
den     = (1+Ka*cafin+Kc*ccfin)*(1+Ka*cafin+Kc*ccfin);
r1      = k1*cafin*cbfin/den;
r2      = k2*ccfin*cbfin/den;
%%
%% adiabatic temperature rise and
%% Ucrit to achieve a zero intial temperature derivative
%%
adrise = -(cafin*delH1+ccfin*delH2)/(Cp*cfin*Mair)
Ucrit  = (r1*delH1 + r2*delH2)/(2/Rt*(Ta-T))

%%
%% global collocation
%%
npts = 40;
[R, A, B, Q] = colloc (npts-2, 'left', 'right');
R = R*Rp;
A = A/Rp;
B = B/(Rp*Rp);
Q = Q*Rp;
Aint = A(2:npts-1,:);
Bint = B(2:npts-1,:);
Rint = R(2:npts-1);

function retval = pellet (x)
  global  npts A Aint Bint Rint  caf cbf ccf ...
          k10 E1 k20 E2 Ka0 Ea Kc0 Ec Da Db Dc T ...
          Ra Rb Rc Da Db Dc kma kmb kmc dcadr dcbdr dccdr
  %%
  %% component A
  %%
  za = x(1:npts);
  zb = x(npts+1:2*npts);
  zc = x(2*npts+1:3*npts);
  ca = exp(za);
  cb = exp(zb);
  cc = exp(zc);

  k1      = k10*exp(-E1/T);
  k2      = k20*exp(-E2/T);
  Ka      = Ka0*exp(-Ea/T);
  Kc      = Kc0*exp(-Ec/T);
  den     = (1+Ka*ca+Kc*cc).^2;
  r1      = k1.*ca.*cb./den;
  r2      = k2*cc.*cb./den;
  Ra      = - r1;
  Rb      = - 1/2*r1  - 9/2*r2;
  Rc      = - r2;
  ip = 1;
  retval(ip,1)    = A(1,:)*za;
  caint = ca(2:npts-1);
  retval(ip+1:ip+npts-2,1) = Bint*za + Aint*za.*(Aint*za + 2./Rint) + ...
      Ra(2:npts-1)./(Da*caint);
  dzadr = A(npts,:)*za;
  retval(ip+npts-1,1) = Da*dzadr - kma*(caf/ca(npts) - 1);
  %%
  %% component B
  %%
  ip = npts+1;
  cbint = cb(2:npts-1);
  retval(ip,1)    = A(1,:)*zb;
  retval(ip+1:ip+npts-2,1) = Bint*zb + Aint*zb.*(Aint*zb + 2./Rint) + ...
      Rb(2:npts-1)./(Db*cbint);
  dzbdr = A(npts,:)*zb;
  retval(ip+npts-1,1) =  Db*dzbdr - kmb*(cbf/cb(npts) - 1);
  %%
  %% component C
  %%
  ip = 2*npts+1;
  ccint = cc(2:npts-1);
  retval(ip,1)    = A(1,:)*zc;
  retval(ip+1:ip+npts-2,1) = Bint*zc + Aint*zc.*(Aint*zc + 2./Rint) + ...
      Rc(2:npts-1)./(Dc*ccint);
  dzcdr = A(npts,:)*zc;
  retval(ip+npts-1,1) = Dc*dzcdr - kmc*(ccf/cc(npts) - 1);
  dcadr = A(npts,:)*ca;
  dcbdr = A(npts,:)*cb;
  dccdr = A(npts,:)*cc;
endfunction
%%
%% find the pellet profile at tube inlet
%%
caf = cafin; cbf=cbfin; ccf=ccfin;
zain = log(1e-9*caf); zaout = log(0.75*caf);
za0 = zain + R/Rp*(zaout-zain);
zbin = log(0.75*cbf); zbout = log(cbf);
zb0 = zbin + R/Rp*(zbout-zbin);
zcin = log(1e-6*ccf); zcout = log(ccf);
zc0 = zcin + R/Rp*(zcout-zcin);

z0 = [za0; zb0; zc0];

[z, fval, info] = fsolve (@pellet, z0)

norm (fval)

plot (R, [z0(1:npts),z(1:npts),z0(npts+1:2*npts),z(npts+1:2*npts), ...
          z0(2*npts+1:3*npts),z(2*npts+1:3*npts)])

_______________________________________________
Bug-octave mailing list
Bug-octave@...
https://www-old.cae.wisc.edu/mailman/listinfo/bug-octave

[Jaroslav Hajek-2]

On Mon, Oct 26, 2009 at 8:20 PM, John W. Eaton <jwe@...> wrote:
> With Octave 3.0.5, fsolve computes a good solution for the following
> problem.  The norm of the function values is approximately 1e-5.  With
> 3.2.3 and Octave from the current hg archive, the norm of the function
> values is about 3500.
>
> jwe
>
>

Er, yes. Believe it or not, in a theoretical sense this is a superior
behavior :)
Here, fsolve reacts to the fact that the problem is very sensitive. Try adding

zerr = z .* (1e-8 * (randn (length (z), 1)));
norm (pellet (z + zerr))

to the script. By introducing a *relative* error of order 1e-8, the
residual jumps up by more than 6 orders of magnitude (I get 5.9322e-06
-> 25.963). The new fsolve automatically scales the variables by
jacobian column norms to get somewhat scaling-independent termination
tests. More precisely, the tests are
residual <= TolFun * n * norm (dg .* x)
norm (dg .* step) <= TolX * norm (dg .* x)

the rhs of the first test is an upper bound on the quantity jacobian *
e, where e is an error of the order TolFun * x. Ergo, TolFun ~ eps
should be close to the best precision achievable (but may be a few
magnitueds off, because the upper bound is not accurate).

By default, the tolerances are set to sqrt(eps), which seems a good
compromise to me. You can achieve what you need by lowering TolFun and
TolX as needed (for example, 1e-16 and 1e-12).

This strategy has the advantage that, for instance, if you scale half
of the variables by 100, you still get the same answer.
The original minpack contained a similar strategy for scaling
variables. However, it is switched off by default when called via
hybrj1, which is done in octave 3.0.x, so 3.0.x's fsolve doesn't have
this invariant property. Further, original MINPACK doesn't use a
termination test on the residual at all, except for a test for an
exact zero. So there doesn't seem to be anything corresponding to
TolFun.

Maybe there could be an option to disable the smart scaling and just
force TolFun and TolX to be hard absolute tolerances?
Or should that be the default? I think the current tests are
mathematically superior, because they take into account the
sensitivity of the problem at hand (which the user may not even be
aware of), but maybe they are too hard to understand...


--
RNDr. Jaroslav Hajek
computing expert & GNU Octave developer
Aeronautical Research and Test Institute (VZLU)
Prague, Czech Republic
url: www.highegg.matfyz.cz

_______________________________________________
Bug-octave mailing list
Bug-octave@...
https://www-old.cae.wisc.edu/mailman/listinfo/bug-octave