Matlab is plotting the graph incorrectly

I have a list of data and some initial parameter estimates (p), and the model i want to use RNp. However, after plotting both (t,RN) and then (t,RNp) on excel using the intial paramter estimates, the curve for (t,RNp) is competely off in Matlab and i cannot understand why.For some reason the values for t are fine but the values for RNp remain constant. I'm fairly new to Matlab so maybe i'm missing something in the coding. And would fmincon or lsqcurvefit be better for the parameter fitting? Thank you for any help :)

[Walter Roberson]

Sorry, what is RNp for this purpose?

[James Thorpe]

It is the reaction rate model for RN with the initial parameter estimates in

[Walter Roberson]

Unfortunately when I copy the image of your code into MATLAB, it just gives me an image of an error message telling me that it can only execute real code, not pictures of code.

Answer 1

[Walter Roberson]

CPA is a vector. It appears on the right hand side of the / expression, so the right hand side is a vector. When you use the / operation with a vector right hand side, you are not invoking division: you are invoking the least-squared fitting operator whose formal name is "matrix right divide".

To get division instead of model fitting, you need to use the ./ operator instead of the / operator.

Note that with your numerator and denominator both being vectors, using the ./ operator would give you a vector result for RNp .

If for some reason you intend matrix-right-divide instead of element-by-element division then at the very least you should add a comment about that. Better would be to re-arrange the terms and use the \ operator: A/B for real matrix A and B is B.'\A.' which MATLAB programmers would immediately recognize as being a fitting instead of a division operator.

[James Thorpe]

Thank you i used the ./ operator and the graph is plotting correctly now.

I have copied and pasted the code so hopefully it will work now, i have tried to fit the model to the data and i thought intiallly that it had worked, the blue curve represents the raw data, the red curve represents the model curve with the intial guess parameters and the green curve should represent the model curve with the fitted parameter. However, the parameter (p) values are not changing and i get the message:

'Local minimum possible. Constraints satisfied.

fmincon stopped because the size of the current step is less than

the value of the step size tolerance and constraints are

satisfied to within the value of the constraint tolerance.'

t =[0, 19, 36, 55, 85, 114, 267];

RN =[0.00001,4.09E-02,4.30E-02,3.56E-02,2.38E-02,1.09E-02,1.17E-03];

CN =[3.44, 2.67, 1.93, 1.26, 0.55, 0.23, 0.05];

CPA =[0, 0.78, 1.52, 2.19, 2.89, 3.2, 3.37];

CSA =[0, 0.0115149336611702, 0.0257909880562478, 0.0386730579526147, 0.0603424603242075, 0.0717329017967482, 0.0796311089653999];

CH =[0.0,5E+00,7E+00,9E+00,12E+00,2E+01,1E+01];

CNH3 =[1, 1, 1, 1, 1, 1, 1];

% p0=[k1 kN kPA kNH3 kH]

p = [8.66E8, 2.515E-6, 3.366E-6, 7.700E-3, 5.033E-6];

RNp = @(p) ((p(1)*p(2)*p(5).*CH.*CN)./((1+p(2).*CN+((p(5).*CH).^0.5)+p(3).*CPA+p(4).*CNH3).^3));

objective = @(p) sum(((RNp(p)-RN)./RN).^2);

%constraints

A = [];

b = [];

Aeq = [];

beq = [];

%bounds

lb = zeros(size(p));

ub =[];

RNp0 = fmincon(objective,p,A,b,Aeq,beq,lb,ub);

disp(['initial objective:' num2str(objective(p))])

disp(['final objective:' num2str(objective(RNp0))])

plot(t,RN,'-bX')

hold on

plot(t,RNp(p),'-ro')

hold on

plot(t,RNp(RNp0),'-go')

disp(['Final p:' num2str((p))])

ylabel('Rate')

xlabel('time')