Non linear multiple curve fitting

Dear all,

I'm trying to do the multiple nonlinear regression of some rheology curves (rheology.png). In particular, viscisty (visc.txt) changes with both temperature (T.txt) and shear rate (shear_rate.txt).

I would find a nonlinear regression of the form given in the attached 'regression_model.txt', where eta is the viscosity, T the temperature and gamma the shear rate.

but I'm having problems... I tried with nonlinfit, but this worked on only if a fix a given value of temperature (so, for a single variable the regression was ok).

How can I perform this regression in a simple way?

I would thanks in advance anyone who will contribute to help me

Best regards,

AP

Answer 1

[Torsten]

Form a big matrix A:

First column: a vector of 1's

Second column: log(gamma) or log10(gamma) (depending on what "lg" means)

Third column: (log(gamma)).^2 or (log10(gamma)).^2 (depending on what "lg" means)

Fourth column: T*log(gamma) or T*log10(gamma) (depending on what "lg" means)

Fifth column: T

Sixth column: T.^2

and a big vector b, consisting of one column with the corresponding values for log(eta) or log10(eta) ( (depending on what lg means)

Then you can solve for the vector v =[ A0; A1; A11; A12; A2; A22] of constants in your model by typing

v = A\b.

[Alessio Pricci]

Hi Torsten,

I tired to implement your suggestion, but it seems that it does not work. Here I report the implementation of the code:

%% Parameters used for the viscosity curves

D1=7.40688e+12;

D2=373.15;

A1=30.62;

A2=51.6;

n=0.2411;

tau=72297.9;

%% shear rate vector (equal for all 3 temperatures)

shear_rate=linspace(1e-1,1e4,1e5);

%% temperature and corresponding viscosity vector

T=210+273.15;

v1=D1*exp(-(A1*(T-373.15))./(A2+(T-373.15)))./((1+((D1*exp(-(A1*(T-373.15))./(A2+(T-373.15)))).*shear_rate./tau)).^(1-n));

T=245+273.15;

v2=D1*exp(-(A1*(T-373.15))./(A2+(T-373.15)))./((1+((D1*exp(-(A1*(T-373.15))./(A2+(T-373.15)))).*shear_rate./tau)).^(1-n));

T=280+273.15;

v3=D1*exp(-(A1*(T-373.15))./(A2+(T-373.15)))./((1+((D1*exp(-(A1*(T-373.15))./(A2+(T-373.15)))).*shear_rate./tau)).^(1-n));

%% Implementation of your suggestion

T1=ones(1,length(shear_rate))*(210+273.15);

T2=ones(1,length(shear_rate))*(245+273.15);

T3=ones(1,length(shear_rate))*(280+273.15);

A1=[ones(1,length(shear_rate)); log(shear_rate); (log(shear_rate)).^2; T1.*log(shear_rate); T1; T1.^2]';

A2=[ones(1,length(shear_rate)); log(shear_rate); (log(shear_rate)).^2; T2.*log(shear_rate); T2; T2.^2]';

A3=[ones(1,length(shear_rate)); log(shear_rate); (log(shear_rate)).^2; T3.*log(shear_rate); T3; T3.^2]';

A=[A1; A2; A3];

b=[log(v1) log(v2) log(v3)]';

v=A\b;

%% Verify at the higher temperature (280[°C])

T=280+273.15;

aaaa=v(1)+v(2)*log(shear_rate)+v(3).*(log(shear_rate).^2)+v(4)*T.*log(shear_rate)+v(5)*T+v(6)*T^2;

aaaa=exp(aaaa);

%% Graphical representation

loglog(shear_rate,aaaa,'b-',shear_rate,v3,'r-');

[Torsten]

I was talking about columns, not rows.

ones(1,length(shear_rate)), e.g., is a row vector, not a column vector.

[Alessio Pricci]

I modified, by transposing the vectors, but the results does not change:

T1=ones(1,length(shear_rate))'*(210+273.15);

T2=ones(1,length(shear_rate))'*(245+273.15);

T3=ones(1,length(shear_rate))'*(280+273.15);

A1=[ones(1,length(shear_rate))', log(shear_rate)', (log(shear_rate)).^2', T1.*log(shear_rate)', T1, T1.^2];

A2=[ones(1,length(shear_rate))', log(shear_rate)', (log(shear_rate)).^2', T2.*log(shear_rate)', T2, T2.^2];

A3=[ones(1,length(shear_rate))', log(shear_rate)', (log(shear_rate)).^2', T3.*log(shear_rate)', T3, T3.^2];

A=[A1; A2; A3];

b=[log(v1) log(v2) log(v3)]';

v=A\b;

T=280+273.15;

aaaa=v(1)+v(2)*log(shear_rate)+v(3).*(log(shear_rate).^2)+v(4)*T.*log(shear_rate)+v(5)*(T)+v(6)*T^2;

aaaa=exp(aaaa);

loglog(shear_rate,aaaa,'b-',shear_rate,v3,'r-');

[Torsten]

D1 = 7.40688e+12;
D2 = 373.15;
A1 = 30.62;
A2 = 51.6;
n = 0.2411;
tau = 72297.9;
T1 = 210+273.15;
T2 = 245+273.15;
T3 = 280+273.15;
shear_rate=(linspace(1e-1,1e4,1e5)).';
V = @(T) D1*exp(-(A1*(T-373.15))./(A2+(T-373.15)))./((1+((D1*exp(-(A1*(T-373.15))./(A2+(T-373.15)))).*shear_rate./tau)).^(1-n));
T = T1;
v1 = V(T1); 
T = T2;
v2 = V(T2); 
T = T3;   
v3 = V(T3); 
A = [ones(3*numel(shear_rate),1),...
    [log10(shear_rate);log10(shear_rate);log10(shear_rate)],...
    [(log10(shear_rate)).^2;(log10(shear_rate)).^2;(log10(shear_rate)).^2],...
    [T1*log10(shear_rate);T2*log10(shear_rate);T3*log10(shear_rate)],...
    [T1*ones(numel(shear_rate),1);T2*ones(numel(shear_rate),1);T3*ones(numel(shear_rate),1)],...
    [T1^2*ones(numel(shear_rate),1);T2^2*ones(numel(shear_rate),1);T3^2*ones(numel(shear_rate),1)]];
b =  [log10(v1);log10(v2);log10(v3)];
v0 = A\b
v = lsqnonlin(@(x)fun(x,A,b),v0)
T = T3;
aaaa = v(1) + v(2)*log10(shear_rate) + v(3)*(log10(shear_rate)).^2 + v(4)*T.*log10(shear_rate) + v(5)*T + v(6)*T^2;
aaaa = 10.^(aaaa);
%% Graphical representation
loglog(shear_rate,aaaa,'b-',shear_rate,v3,'r-');

function res = fun(x,A,b)
  res = 10.^(A*x) - 10.^b;
end

[Torsten]

Be careful with the log-expressions:

https://en.wiktionary.org/wiki/lg