fitting is not working

%% Luz-Meiboom fit of R2 vs pH
clear; clc; close all;
%% Data
pH = [1.88, 2.45, 3.39, 3.87, 4.53, 4.68, 5.12, 5.25, 5.55, 5.82, ...
    6.01, 6.27, 6.66, 6.98, 7.75, 8.25, 8.96, 9.23];
R2 = [0.01137, 0.00973, 0.00803, 0.01824, 0.0547, 0.0712, 0.10246, ...
    0.12065, 0.12992, 0.22537, 0.18405, 0.05077, 0.0203, 0.01554, ...
    0.02122, 0.0371, 0.03736, 0.00611];
pH = pH(:);
R2 = R2(:);
%% Luz-Meiboom model function
% params(1) = R2m
% params(2) = tau
luzMeiboomFun = @(params, x) params(1) + params(2).* ...
    (1 - (2.*params(2)./x).*tanh(x./(2.*params(2))));
%% Initial parameter guesses
R2m_guess  = min(R2);               
tau_guess  = 0.963;                     
p0 = [R2m_guess, tau_guess];
%% Bounds (lower, upper)
lb = [0,    0,    1e-4];
ub = [1,    1000,  500];
%% Fit using lsqcurvefit (Optimization Toolbox)
opts = optimoptions('lsqcurvefit', 'Display', 'iter');
[params_fit, resnorm, residual, exitflag] = lsqcurvefit(luzMeiboomFun, p0, pH, R2, lb, ub, opts);

Warning: Length of lower bounds is > length(x); ignoring extra bounds.

Warning: Length of upper bounds is > length(x); ignoring extra bounds.

                                            Norm of      First-order 
 Iteration  Func-count      Resnorm            step       optimality
     0          3           5.69629                             3.32
     1          6           1.34129        0.587386             1.28      
     2          9           0.32902        0.306607            0.356      
     3         12           0.11218        0.185613           0.0897      
     4         15         0.0749423        0.106142           0.0195      
     5         18         0.0712123       0.0460272          0.00271      
     6         21         0.0711088      0.00972236         0.000108      
     7         24         0.0711056       0.0128557           0.0159      
     8         27         0.0711049      0.00688321          0.00627      

Local minimum possible.

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

R2m_fit  = params_fit(1);
tau_fit  = params_fit(2);
fprintf('\n--- Fit results ---\n');

--- Fit results ---

fprintf('R2m  = %.5f\n', R2m_fit);

R2m  = 0.00000

fprintf('tau  = %.5f\n', tau_fit);

tau  = 0.06460

fprintf('Residual norm = %.5f\n', resnorm);

Residual norm = 0.07110

%% Generate smooth fitted curve for plotting
pH_fine = linspace(min(pH), max(pH), 300);
R2_fit_curve = luzMeiboomFun(params_fit, pH_fine);
%% Plot: pH vs R2 
figure('Color','w');
plot(pH, R2, 'ko', 'MarkerFaceColor','r', 'MarkerSize',7, 'DisplayName','Data');
hold on;
plot(pH_fine, R2_fit_curve, 'b-', 'LineWidth', 1.8, 'DisplayName','Luz-Meiboom fit');
xlabel('pH', 'FontSize', 12);
ylabel('R_2 (s^{-1})', 'FontSize', 12);
title('R_2 vs pH — Luz-Meiboom Fit', 'FontSize', 13);
legend('Location','best');
grid on;
box on;
hold off;

[Walter Roberson]

Your upper and lower bounds are three elements each, but you are only searching over 2 parameters.

Answer 1

[Torsten]

As far as I can see, the function f(x) = p*(1-2*p./x.*tanh(x/(2*p))) is monotonically increasing (the first parameter params(1) is not important in the analysis since it's only a translation of the curve in the upward or downward direction). Thus you will never get the bellshape of your experimental data for a fit result. Try plotting the function from above for different values of p, and you will see what I mean.

Thus in short: the model function is not suited to reflect your experimental data.

x = 1:14;
p = [0.003;0.03;0.3;3;30;300];
f = p.*(1-2*p./x.*tanh(x./(2*p)));
plot(x,f)

[John D'Errico]

It appears you (both) may be using the wrong model form. That is, when I look online for that model, I see a different form.

https://www.nmr-relax.com/manual/The_LM63_3_site_fast_exchange_CPMG_model.html

In there (aswell in all other places I saw that model), the form appears to have an extra 1/x term in it. Something more like this (again, I've dropped extraneous constants so we can see the fundamental shape.)

f = @(x,p) 1./x.*(1-2*p./x.*tanh(x./(2*p)));
figure
fplot(@(x) f(x,1),[0,50],'r')
hold on
fplot(@(x) f(x,2),[0,50],'g')
fplot(@(x) f(x,3),[0,50],'b')

Anyway, this now has a skewed bell shape to it. At the same time, the shape is still not terribly consistent with the data provided, so I doubt it will help greatly.

[Ehtisham]

Can you rewrite the code with this model

[Torsten]

Replace

luzMeiboomFun = @(params, x) params(1) + params(2).* ...
    (1 - (2.*params(2)./x).*tanh(x./(2.*params(2))));

by the modified model function.

[Ehtisham]

% Luz-Meiboom fit of R2 vs pH
clear; clc; close all;
%% Data
pH = [1.88, 2.45, 3.39, 3.87, 4.53, 4.68, 5.12, 5.25, 5.55, 5.82, ...
      6.01, 6.27, 6.66, 6.98, 7.75, 8.25, 8.96, 9.23];
R2 = [0.01137, 0.00973, 0.00803, 0.01824, 0.0547, 0.0712, 0.10246, ...
      0.12065, 0.12992, 0.22537, 0.18405, 0.05077, 0.0203, 0.01554, ...
      0.02122, 0.0371, 0.03736, 0.00611];
pH = pH(:);
R2 = R2(:);
%% Luz-Meiboom model function
% params(1) = R2m
% params(2) = tau
luzMeiboomFun = @(params, x) params(1) + params(2).* ...
    (1 - (2.*params(2)./x).*tanh(x./(2.*params(2))));
%% Initial parameter guesses
R2m_guess  = min(R2);                 
tau_guess  = 0.963;                    
p0 = [R2m_guess,  tau_guess];
%% Bounds (lower, upper) 
lb = [0,    0,    1e-4];
ub = [1,    100,  50];
%% Fit using lsqcurvefit 
opts = optimoptions('lsqcurvefit', 'Display', 'iter');
[params_fit, resnorm, residual, exitflag] = lsqcurvefit(luzMeiboomFun, p0, pH, R2, lb, ub, opts);

Warning: Length of lower bounds is > length(x); ignoring extra bounds.

Warning: Length of upper bounds is > length(x); ignoring extra bounds.

                                            Norm of      First-order 
 Iteration  Func-count      Resnorm            step       optimality
     0          3           5.69629                             3.32
     1          6           1.34129        0.587386             1.28      
     2          9           0.32902        0.306607            0.356      
     3         12           0.11218        0.185613           0.0897      
     4         15         0.0749423        0.106142           0.0195      
     5         18         0.0712123       0.0460272          0.00271      
     6         21         0.0711088      0.00972236         0.000108      
     7         24         0.0711056       0.0128557          0.00159      
     8         27         0.0711049      0.00698329         0.000613      

Local minimum possible.

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

R2m_fit  = params_fit(1);
tau_fit  = params_fit(2);
fprintf('\n--- Fit results ---\n');

--- Fit results ---

fprintf('R2m  = %.5f\n', R2m_fit);

R2m  = 0.00000

fprintf('tau  = %.5f\n', tau_fit);

tau  = 0.06461

fprintf('Residual norm = %.5f\n', resnorm);

Residual norm = 0.07110

%% Generate smooth fitted curve for plotting
pH_fine = linspace(min(pH), max(pH), 300);
R2_fit_curve = luzMeiboomFun(params_fit, pH_fine);
%% Plot: pH vs R2 
figure('Color','w');
plot(pH, R2, 'ko', 'MarkerFaceColor','r', 'MarkerSize',7, 'DisplayName','Data');
hold on;
plot(pH_fine, R2_fit_curve, 'b-', 'LineWidth', 1.8, 'DisplayName','Luz-Meiboom fit');
xlabel('pH', 'FontSize', 12);
ylabel('R_2 (s^{-1})', 'FontSize', 12);
title('R_2 vs pH — Luz-Meiboom Fit', 'FontSize', 13);
legend('Location','best');
grid on;
box on;
hold off;

i replace that but giving the same plot

[Torsten]

I don't see changes in

luzMeiboomFun = @(params, x) params(1) + params(2).* ...
    (1 - (2.*params(2)./x).*tanh(x./(2.*params(2))));

compared to your previous code.

And you still use three-element vectors for lower and upper bounds on the parameters instead of two-element vectors.

[Ehtisham]

See i change the funcztion but still getting the same result

% Luz-Meiboom fit of R2 vs pH
clear; clc; close all;
%% Data
pH = [1.88, 2.45, 3.39, 3.87, 4.53, 4.68, 5.12, 5.25, 5.55, 5.82, ...
      6.01, 6.27, 6.66, 6.98, 7.75, 8.25, 8.96, 9.23];
R2 = [0.01137, 0.00973, 0.00803, 0.01824, 0.0547, 0.0712, 0.10246, ...
      0.12065, 0.12992, 0.22537, 0.18405, 0.05077, 0.0203, 0.01554, ...
      0.02122, 0.0371, 0.03736, 0.00611];
pH = pH(:);
R2 = R2(:);
%% Luz-Meiboom model function
% params(1) = R2m
% params(2) = tau
luzMeiboomFun = @(params, x) params(1) + params(2).* ...
    (1 - (2.*params(2)./x).*tanh(x./(2.*params(2))));
%% Initial parameter guesses
R2m_guess  = min(R2);                 
tau_guess  = 1;                    
p0 = [R2m_guess,  tau_guess];
%% Bounds (lower, upper) 
lb = [0,   1e-4];
ub = [1,   50];
%% Fit using lsqcurvefit 
opts = optimoptions('lsqcurvefit', 'Display', 'iter');
[params_fit, resnorm, residual, exitflag] = lsqcurvefit(luzMeiboomFun, p0, pH, R2, lb, ub, opts);

                                            Norm of      First-order 
 Iteration  Func-count      Resnorm            step       optimality
     0          3           5.94767                             3.34
     1          6           1.39853        0.604756             1.32      
     2          9          0.342046         0.31134             0.37      
     3         12          0.114716        0.188396           0.0935      
     4         15         0.0752786        0.108196           0.0205      
     5         18         0.0712276       0.0475836          0.00291      
     6         21         0.0711087       0.0103039         0.000123      
     7         24         0.0711056       0.0122266         0.000597      
     8         27         0.0711049      0.00718406          0.00032      

Local minimum possible.

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

R2m_fit  = params_fit(1);
tau_fit  = params_fit(2);
fprintf('\n--- Fit results ---\n');

--- Fit results ---

fprintf('R2m  = %.5f\n', R2m_fit);

R2m  = 0.00000

fprintf('tau  = %.5f\n', tau_fit);

tau  = 0.06461

fprintf('Residual norm = %.5f\n', resnorm);

Residual norm = 0.07110

%% Generate smooth fitted curve for plotting
pH_fine = linspace(min(pH), max(pH), 300);
R2_fit_curve = luzMeiboomFun(params_fit, pH_fine);
%% Plot: pH vs R2 
figure('Color','w');
plot(pH, R2, 'ko', 'MarkerFaceColor','r', 'MarkerSize',7, 'DisplayName','Data');
hold on;
plot(pH_fine, R2_fit_curve, 'b-', 'LineWidth', 1.8, 'DisplayName','Luz-Meiboom fit');
xlabel('pH', 'FontSize', 12);
ylabel('R_2 (s^{-1})', 'FontSize', 12);
title('R_2 vs pH — Luz-Meiboom Fit', 'FontSize', 13);
legend('Location','best');
grid on;
box on;
hold off;

[Walter Roberson]

Your code

luzMeiboomFun = @(params, x) params(1) + params(2).* ...
    (1 - (2.*params(2)./x).*tanh(x./(2.*params(2))));

does not contain a 1./x term.

[Ehtisham]

@Torsten mention just to chnage this function like this luzMeiboomFun = @(params, x) params(1) + params(2).* ...

(1 - (2.*params(2)./x).*tanh(x./(2.*params(2))));