optimize fsolve for diffusion and reaction in a spherical particle

Hello to everyone,

I am trying to solve a system of nonlinear equations describing diffusion and second order reaction in a spherical catalyst particle. The equations are formed by making the problem dimensionless and using the finite differences method. The problem is that even though the resulting plot seems about right I always get an exit flag of 2 or 3. I really don't know a lot about matlab and the information on optimization I found so far seem overwelming. I used some of the suggested values but no lack so far. Are the equations not defined properly? Please have a look at the code below:

function catal_partF

clc

clear all

global theta h Thiele Biot c N

N = 1001;

h = 1/(N-1);                 %step size

theta = 0:h:1;               %dimensionless legth

p = length(theta);

Biot = [100 10 1 0.1 0.01];

for i = 1:5
    Bi = Biot(i);

      q = 1;
      for Thiele = logspace(-1,2,10)
          c0 = linspace(0.5,1,p);           %initial values for the concentration           
          options = optimset('Display','iter','Tolfun',1e-10,'TolX',1e-10,'ScaleProblem','Jacobian');
          [c,fval,exitflag,output] = fsolve(@sol,c0)
          eff_num(q,i) = trapz(theta,3.*(c.^2).*theta.^2) %numerical value of the effectiveness factor
          q = q+1;
      end
  end

Thiele = logspace(-1,2,10);
figure(5)
loglog(Thiele,eff_num)
title('Effectiveness factor versus the Thiele modulus','FontWeight','bold')
xlabel('Thiele modulus','FontWeight','bold'), ylabel('effectiveness factor','FontWeight','bold')
legend('Bi=100','Bi=10','Bi=1','Bi=0.1','Bi=0.01')

function F = sol(c)

F(1)=(-3*c(1)./(2*h))+(4*c(2)./(2*h))-(c(3)./(2*h));
F(2:N-1)=(1/h^2-1./(theta(2:N-1)*h)).*c(1:N-2)-(2/h^2)*c(2:N-1)-Thiele^2*c(2:N-1).^2+(1/h^2+1./(theta(2:N-1)*h)).*c(3:N);
F(N)=(3/(2*h)+Bi).*c(N)-(4/(2*h)).*c(N-1)+c(N-2)./(2*h)-Bi;

end

end

Thank you in advance!

Answer 1

[Alan Weiss]

It is difficult to know what you are looking for. Your code runs, so it is at least sensible to MATLAB. And, generally speaking, an exit flag of 2 or 3 is not too bad, fsolve thinks it found a good answer in that it reduced the value of the vector function to a small enough value. It did not give an exit flag of one because it could not reduce the estimated derivative enough.

I am a little surprized that you use trapz instead of one of the more efficient integration methods, such as ode45.

And it seems to me that you are integrating your answer, but you don't start fsolve from the previous answer. You might try doing just that; presumably, the solution that fsolve gets doesn't change much froom integration point to integration point, so it might be good to start with the last solution, not the same point c0 every time.

Also, for the sake of efficiency, do not call optimset inside a loop. Set your options once, at the beginning of your calculations, and reuse them.

Good luck,

Alan Weiss

MATLAB mathematical toolbox documentation

[Nikoleta]

Hello and thank you for your answer,

You are of course right about the optimset command.

I use the trapz command instead of for example ode45 only because I am asked to, you see this is supposed to be part of an assignment.

Actually the results here are not bad but I know that something is wrong, either with the setup of the functions (but I don't think it is that), or the layout of the code and the initial estimations (c0).

Here is why; In this case you can see that I have a fixed number of grid points, N, for the discetization and a varied Biot number. In another case I use more or les the same code but the N varies and the Biot number is fixed. There the problem is more obvious because the result (c, the concentration) is disturbingly low and doesn't even come close to the value of 1 for theta=1 as it is supposed to.

I must optimize my code and your point on using the last result each time as an estimation for the next iteration is correct. To be honest I thought about it but I have a bit of a struggle as to how to integrate this in my code..