problem with fsolve - unpredictable results with long data vectors

Hello all. I am using fsolve to solve "fraction bound" as part of larger function that analyzes calorimetry data, though that isn't important. I have a function that describes the binding (known as the partition function or binding polynomial). It is simple to write, but finding exact (analytical) solutions can be hard. Hence, I intend on using fsolve. I am starting with the simplest case, which does have an easy solution, but later ones will be more complex.

My problem is that the output depends on the size of my data. For short data sets (limited number of injections" the output "fraction bound" is continuous and looks like it should. However, if I use long data sets (lots of injections), I get strange behavior; there will be a few values that are not continuous and basically wrong. It is probably some sort of minimization issue. Whilst playing around, I found that tinkering with the finite-difference settings helped. So, I wrote my own Jacobian. Surprisingly this didn't help.

I am very interested in anyone's take as to what is happening. I am inserting two files. One is a function and the other a script that calls the function.

function [ X ] = testfun( X0, i, Mt, Xt, n, Ka )
options=optimset(...
    'Algorithm','trust-region-dogleg',...
    'Display','iter',...
    'DerivativeCheck','on',...
    'Jacobian','on',...
    'Diagnostics','on',...
    'DiffMaxChange',0.01,...
    'FinDiffType','central',...
    'FunValCheck','on',...
    'MaxIter',1e3,...
    'TolFun',1e-12,...
    'TolX',1e-9,...
    'TypicalX',0.1.*ones(1,length(i))');

      function [BP, J] = nested_BP(X)
          Lf = Xt - X.*Mt;
          Z = (1 + Ka.*Lf).^n; 
          dZf = Ka.*n.*(1 + Ka.*Lf).^(n-1); 
          BP = (dZf.*Lf)./Z - X; 
          dLf = -Mt;
          dZx = dZf .* dLf; 
          ddZfx = (n-1).*n.*Ka.^2.*(1 + Ka.*Lf).^(n-2).*dLf; 
          % d(BP)/dX = d(dZf.*Lf)./Z + (dZf.*Lf)d(Z^-1) - 1
          dBPx = (ddZfx.*Lf + dZf.*dLf)./Z - (dZf.*Lf).*(Z.^-2).*dZx - 1;
          NN = length(i);
          J = sparse(1:NN,1:NN,dBPx,NN,NN);
      end
  X = fsolve(@nested_BP, X0, options);
end

Now the script:

clear;
Mt0 = 10e-6;
V0 = 1e-3;
Vi = 10e-6;
Xi = 0.1e-3;
n = 4;
Kd = 1e-6;
Ka = 1/Kd;
i = 1:200;
Mti = Mt0 .* ( 1 - (Vi/V0) ).^i;
Xti = Xi .* ( 1 - ( 1 - (Vi/V0) ).^i);
X0 = Kd.*n.*Xti./(Mti); 
X = testfun(X0,i,Mti,Xti,n,Ka);
plot(i,X);

As written, it works just fine. You can see how the output is a smoothly increasing function with an asymptote approaching 4.

Now, change from 200 to 400 injections:

i=1:400

And try it again. See all those crazy values in the first 50 or so injections? I am trying to understand where they come from. I am guessing that the problem lies in the all the extra values that approach the asymptote. But, I can't predict in general, what sort of data I am going to be analyzing. And this function is going to part of a larger non-linear least squares fitting function. Having this sort of thing happening unpredictably during the fitting may lead to disaster (or not).

Answer 1

[Steve Grikschat]

I'm not sure what exactly is going on here. Fsolve seems to be finding a solution in both cases, correct? Since the system is square, there is only one solution. From what I can tell, then the instability must be a property of the model (or its implementation), but I don't know enough about it to pinpoint the cause.

Just a guess, but maybe the problem related to the scaling introduced when you increase the number of and magnitude of the values in "i". In your setup, you take values to the power of i:

Mti = Mt0 .* ( 1 - (Vi/V0) ).^i;
Xti = Xi .* ( 1 - ( 1 - (Vi/V0) ).^i);

The larger the value of i, the smaller some of these values get. Could that be it?

+Steve