How to improve precision of lsqnonlin analysis

I used lsqnonlin annlysis to fit a model to two sets of data. The results showed that the differences between original data and fitted line are rather huge (see in Figure1). How to improve this analysis to obtain better results? Thanks! Below are my codes.

function Result = myfunc1(a,p,G,m,fraction)
totalconcentration = 0.0000254;
n = 2;
K = exp(-((G + m .* fraction)./(8.314*298)));
opts = optimoptions('fsolve','Algorithm', 'levenberg-marquardt','FunctionTolerance',1.0000e-12);
monomerconcentration = fsolve(@denaturationfun,0.999,opts);
function y = denaturationfun(x)
    y = K.*totalconcentration - a.^(-1).*(((a.*x).^(n+1)).*(n.*a.*x-n-1)./((a.*x-1).^2)+a.*x./((a.*x-1).^2))+a.^(n-1).*((x^(n+1)).*(n.*x-n-1)./((x-1).^2));
end
Result = p * (1- (monomerconcentration ./ ( K * totalconcentration )));
end


clc;
clear;
close all;
fractiondata = [0.0431	0.0478	0.0525	0.0571	0.0617	0.0662	0.0707	0.0751	0.0795	0.0839	0.0882	0.0925	0.0967	0.1009	0.1051	0.1092	0.1133	0.1174	0.1253];
DegreeofAggdata = [0.86089	0.90051	0.84268	0.9543	0.98855	0.98538	1	0.98493	0.91339	0.92209	0.85817	0.78529	0.64172	0.45712	0.24855	0.11291	0.00812	0	0.0169];
fun = @(x)DegreeofAggdata - myfunc1(x(1),x(2),x(3),x(4),fractiondata);
x0 = [0.001,1,-50000,100000];
lb = [0,0.7,-150000,20000];
ub = [1,1.3,-5000,150000];
options = optimoptions('lsqnonlin','Algorithm','levenberg-marquardt');
R = lsqnonlin(fun,x0,lb,ub);      %% fit model to exp data
A = 0.03:0.001:0.15;
D = zeros(length(A),1);
for i = 1:length(A)
    y = myfunc1(R(1),R(2),R(3),R(4),A(i));
    D(i)= y;
end   
figure(1);
plot(A,D,'r-');      %% draw lines with fitted parameters
hold on;
plot (fractiondata,DegreeofAggdata,'-go');

[Torsten]

monomerconcentration = fsolve(@denaturationfun,0.999,opts);

You only specify one initial value (0.999) for "fsolve". That means that you will get back one value for "monomerconcentration". In my understanding, you want to get back values for "monomerconcentration" for each value in the array K (thus as many values as the array K has elements). Or am I wrong ?

[Haotian]

Yes, K is derived from the 'fraction' data array and i want to get back different 'monomerconcentration' of each 'fraction' data in the array. So that means my function get only one 'monomerconcentration' in the whole process? Do i need to add loop or something in the function file to correct this?

[Shashi Kiran]

@Haotian ,

This might help.

function Result = myfunc1(a, p, G, m, fraction)
    totalconcentration = 0.0000254;
    n = 2;
    opts = optimoptions('fsolve', 'Algorithm', 'levenberg-marquardt', 'FunctionTolerance', 1.0000e-12, 'Display', 'off');
    monomerconcentration = zeros(size(fraction)); 

    for i = 1:length(fraction)
        K = exp(-((G + m * fraction(i)) / (8.314 * 298)));
        monomerconcentration(i) = fsolve(@(x) denaturationfun(x, a, n, K, totalconcentration), 0.999, opts);
    end

    Result = p * (1 - (monomerconcentration ./ (exp(-((G + m .* fraction) ./ (8.314 * 298))) * totalconcentration)));
end

function y = denaturationfun(x, a, n, K, totalconcentration)
    y = K * totalconcentration - a^(-1) * (((a * x)^(n + 1)) * (n * a * x - n - 1) / ((a * x - 1)^2) + a * x / ((a * x - 1)^2)) + a^(n - 1) * ((x^(n + 1)) * (n * x - n - 1) / ((x - 1)^2));
end

[Matt J]

@Haotian The inner problem that you are solving with fsolve is a search for a root of a rational function of x. As you know, a rational function will in general have multiple discrete roots. Your code doesn't appear to have a process for deciding which of these roots to select.

[Torsten]

The solution of the inner problem using "fsolve" does not seem to work in each case. You should try to make it more reliable. And - as @Matt J mentionned - there might be multiple solutions, and "fsolve" could pick the wrong one.

clc;
clear;
close all;
fractiondata = [0.0431	0.0478	0.0525	0.0571	0.0617	0.0662	0.0707	0.0751	0.0795	0.0839	0.0882	0.0925	0.0967	0.1009	0.1051	0.1092	0.1133	0.1174	0.1253];
DegreeofAggdata = [0.86089	0.90051	0.84268	0.9543	0.98855	0.98538	1	0.98493	0.91339	0.92209	0.85817	0.78529	0.64172	0.45712	0.24855	0.11291	0.00812	0	0.0169];
fun = @(x)DegreeofAggdata - myfunc1(x(1),x(2),x(3),x(4),fractiondata);
x0 = [0.001,1,-50000,100000];
lb = [0,0.7,-150000,20000];
ub = [1,1.3,-5000,150000];
options = optimoptions('lsqnonlin','Algorithm','levenberg-marquardt');
format long
R = lsqnonlin(fun,x0,lb,ub)      %% fit model to exp data

Local minimum possible.

lsqnonlin stopped because the final change in the sum of squares relative to 
its initial value is less than the value of the function tolerance.

R = 1×4

1.0e+05 * 0.000000000005208 0.000010157108440 -0.429540795547637 1.496768116265411

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

A = 0.03:0.001:0.15;
D = zeros(length(A),1);
for i = 1:length(A)
    y = myfunc1(R(1),R(2),R(3),R(4),A(i));
    D(i)= y;
end   
figure(1);
plot(A,D,'r-');      %% draw lines with fitted parameters
hold on;
plot (fractiondata,DegreeofAggdata,'-go');

function Result = myfunc1(a,p,G,m,fraction)
  totalconcentration = 0.0000254;
  n = 2;
  K = exp(-((G + m .* fraction)./(8.314*298)));
  opts = optimoptions('fsolve','Algorithm', 'levenberg-marquardt','Display','off','FunctionTolerance',1.0000e-12);
  monomerconcentration = fsolve(@denaturationfun,0.999*ones(1,numel(K)),opts);
  function y = denaturationfun(x)
      y = K.*totalconcentration - a.^(-1).*(((a.*x).^(n+1)).*(n.*a.*x-n-1)./((a.*x-1).^2)+a.*x./((a.*x-1).^2))+a.^(n-1).*((x.^(n+1)).*(n.*x-n-1)./((x-1).^2));
  end
  Result = p * (1- (monomerconcentration ./ ( K * totalconcentration )));
end