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;
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:
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).
