Hi, I need to determine k and n for each unique temperature from the attached text file.

I have a power function: r = kc^n and I need to find k and n for each unique temperature in the attached 'reactions.txt' text file.

I would like to linearise the data in order to plot the fitted curves for each temperature, so I have attached a linear regression function file that I have written.

The code below is what I have so far but I am struggling to linearise the power function.

Thanks

Answer 1

[Walter Roberson]

r = k*c^n so

log(r) = log(k) + n *log(c)

so per temperature solve the matrix

s=[log(c(:)), ones(numel(c))] \ log(r(:))

and then

n = s(1)
k = exp(s(2))

This is technically a nonlinear regression rather than a linear regression. If it is important to use a linear regression then the values calculated should make very good starting points.

[Cate may]

I am trying to do this with linear regression with the following linreg function file that I've written, however I'm getting k as a complex number which doesn't seem right. My code and function file is below:

function [a0,a1,r2] = linreg(x,y)
% INPUTS:
%  - x: linear independent data set 
%  - y: linear dependent data set 
% OUTPUT:
%  - a0: constant in y=a1*x + a0 
%  - a1: gradient in y=a1*x + a0
%  - r2: coefficient of determination

% Getting best regression coefficients
n = length(x);
sx = sum(x);
sy = sum(y);
sx2 = sum(x.^2); %sx squared
sxy = sum(x.*y); %sum of x times y
a1 = (n*sxy-sx*sy)/(n*sx2-sx^2); 
a0 = mean(y) - a1*mean(x); 

% Getting r2 value
st = sum((y-mean(y)).^2); 
sr = sum((y-a0-a1*x).^2);
r2 = (st-sr)/st;

%m-file is below:
r_data = importdata('reactions.txt');
reaction_data = r_data.data;
t = reaction_data(:,1); %temperatures in vector
c = reaction_data(:,2); % concentration
r = reaction_data(:,3); % reaction rate
% Reaction rate equation
% %k = reaction constant
% %c = concentration
% %n = reaction order. 

 temp = unique(t);
for a = 1:length(temp)
    L= t==temp(a) ;
    if temp(a)==323
        x = log(c);
        y = log(r);
[a0,a1,r2] = linreg(x,y)
    end 
end 
k = log(a0);
n = a1;

[Walter Roberson]

K>> [min(x),max(x)]
ans =
          -1.6094379124341         0.587786664902119

K>> [min(y),max(y)]
ans =
         -3.96859335691654          1.09296302826941
    

a power function: r = kc^n and I need to find k and n

You are processing with n non-integer. A^B when A is non-integer is defined as exp(B*log(A)). log() of a negative number A comes out as 1i*pi + log(-A) . Multiply that by B and you get B*1i*pi + B*log(-A) . exp() of that will be real-valued only in the case where B is an integer. Is your power to be discovered, n, certain to be an integer? No, not at all!! So now the question becomes whether the c is certain to be positive. But as you look at both x and y you can see that they both involve negative quantities and positive quantities, so no matter whether you c is x or your c is y, you are going to be raising a negative quantity to a floating point n.

Is there any hope of rescuing the model? Hypothesize that there is an integer n for r = k*c^n .

Case 1: n is even. Then when c is negative or c is positive, c^n would be positive, and sign(r ) would be sign(k) .

Is c = y and r = x ?

         [x(1),y(1); x(4),y(4)]

ans =

-0.510825623765991 -0.958676448370593

0.182321556793955 -0.127719741602371

negative y to an even integer would be positive, so for x(1) to be negative, k would have to be negative. But if k is negative then with negative y(4) to an even integer would be positive, multiply by negative k and you would get negative for x(4). Therefore, either n is not an even integer or else c = y and r = x is not true.

Is c = x and r = y ?

         [x(1),y(1); x(7),y(7)]

ans =

-0.510825623765991 -0.958676448370593

0.587786664902119 1.09296302826941

negative x(1) to an even integer would be positive, so for y(1) to be negative, k would have to be negative. But if k is negative then with positive x(7) to an even integer being positive, multiply by negative k, and you would get negative for x(4). Therefore, either n is not an even integer or else c = x and r = y is not true.

As we have ruled out c = x and r = y and c = y and r = x for even integer n, we have ruled out even integer n

Case 2: n is odd. Then when c is negative, c^n would be negative, and when c is postive, c^n would be positive, and k*c^n would have sign equal to sign(k)*sign(c ), so sign(x) = sign(y) always for positive k, and sign(x) = -sign(y) always for negative k . Now look back at the displays for indices for 1 and 4, above, and see that for x(1),y(1) that the signs agree (implying positive k), but for x(4),y(4) that the signs disagree (implying negative k). Therefore, n is not an odd integer.

We have now ruled out n being an even integer and have ruled out n being an odd integer.

Therefore your model cannot be fit with integer n, and to have a hope of getting your fitting to not produce a complex k, you would would have to carefully define what it means for c^n when c is negative and n is not an integer. For example if trial n is 1/2 then how do you want to define (-1)^(1/2) in a way that produces a real-valued result ?