Help: another question on estimating reaction rate constants with mathcad

[FanCG]

Hi, everyone

i m a greenhand on mathcad, and I want to solve some of my old problems using mathcad, which have been solved by matlab and maple .

yesterday I encountered a chemical kinetics problem, which I can easily handle through matlab nonlinear toolbox. I reviewed some similar threads in our forum, and cannot apply the method (Minerr, Minimize,etc) into my worksheet.

Maybe I need some help (see the attached mathcad14 or 11 file where I describe my question)

[jeanGiraud]

If you don't get a quick answer, "save as" Mathcad version 11or lower

Just the data set, nothing else.

jmG

[AlanStevens · Accepted solution]

See attached, saved in M14/15 format as M11 doesn't find a solution.

You have k1 as the rate constant for "B". I've assumed this should be k2 - you'll need to change this if I'm wrong.

There is probably a slicker solution but it's a start!

Alan

NB I've edited this file as I realised the solution I posted earlier was wrong.

[FanCG]

Alan , thank you so much !!!!!!

you found a very obvious mistake in my model postulates (it should be k2), which I HADN'T FOUND AFTER TIMES OF CHECK LAST NIGHT. thank you, thank you !!! and you update your reply so fast that my reply cannot keep up with you. while i don't think there are any serious mistakes in your earlier solutions

the "Coeff = 22.4*10^6" in my question is just a coefficient used to convert the units of my data from ' ppmv ' into 'mol/L' at standard condition (because I changed solids volume, mass and reactor volume in my experiment). maybe it is not convient for numerical calculaion. So, in your old worksheet, the unit of raw data is actually converted into kmol/L .

a (ppmv) => a (10^ -6) / 22.4(L/mol) = [a / (22.4 * 10^6) ] (mol/L) = [a / (22.4 * 10^3) ] (kmol/L)

btw, parameter 'm' in my question is a explicit function of bulk density, and has a unit of g/L (or kg/m3, actually they are the 'same') according to the unit of gas concentration, so the rate constant unit you got is

e.g. k2: m3_gas/(kg_solid * day) or L/(g * day) ,

but for k1, a little more complicated, will be [kmol^0.5 * L^0.5 * g^-1 *day^-1]

your results for k are

k1 = 6.211 E-05 [kmol^0.5 * L^0.5 * g^-1 *day^-1]
k2 = 5.993 E-05 L/(g day)
k3 = 4.080 E-04 L/(g day)

and your new results are

k1 = 1.878 E-06 [mol^0.5 * L^0.5 * g^-1 *day^-1]
k2 = 6.323 E-05 L/(g day)
k3 = 4.081 E-04 L/(g day)

and my old result for k (in matlab, using lsqnonlin) is

k1 = 1.24 E-06 [mol^0.5 * L^0.5 * g^-1 *day^-1]
k2 = 3.84 E-05 L/(g day)
k3 = 2.66 E-04 L/(g day)

and I think the earlier results make sense.

and the difference between mathcad and matlab is acceptable compared with the difference between the results from different temperatures. I will try some other operating conditions in mathcad......

[AlanStevens]

I thought my first set of results were invalid because of the square roots of c. It meant the equations didn't all scale the same way. What units are assumed for each parameter in the equations? Any scaling factor would cancel out of the rate equation for "B", but would not for the other two, hence skewing the fits.

Are you absolutely sure of the square roots? If you replace the square roots of c by just c you get a much better fit!!

Alan

[FanCG]

my postulates on the reaction scheme are that :

1) dA/dt = k1 *m* sqrt(C)

2) dB/dt = k2 * m* C

dC/dt = - k3 *m* C (Note that dC/dt here is only written for an assumed adsorption process(C --> C*|solids), it is a rough estimation, i.e., dC/dt here only accounts for a assumed physical adsorption, the physical significant of k3 may change after comparing all the results from different temperatures, and this equation will be incorporated into the following equation)

3) -dC/dt = k1 *m* sqrt(C) + k2 * m *C + k3 * m * C

for reaction 1), I guess it is a half order reaction from other experimental data, and actually, it should be even less than half an order process

in addition, the reactions listed above are not written according to real, or pure chemical equations we learn in textbooks, which can be easily balanced with various method. The solid is a very complex mixture, and roughly viewed as a whole to participate in the reactions.

I noted that you use dimensionless sum of square, it is a great help for me, I directly use Σ(f-fexp)^2 in my earlier calculation as I did in matlab, you let me recall a good manner.

and I will study your new worksheet more carefully after work, thanks a lot for your warm hearted help.

[AlanStevens]

1. Fan CG

"my postulates on the reaction scheme are that :

1) dA/dt = k1 *m* sqrt(C)"

I wonder about the validity of this postulate given that the fits to the empirical data seem so much better with sqrt(C) replaced by C. Is there a theoretical reason to expect something close to the square root behaviour?

2. jean

It would be interesting to see your solution method (I didn't actually check if my latest solution worked in M11, nor did I attempt to tweak the first one to see if I could get it to find one). It's probably much more straightforward than mine (which seemed rather cumbersome).

Alan

[AlanStevens]

Just to reinforce the question over the validity of the sqrt(C) assumption (which jean's method might show to be a non-question] here are the comparisons with sqrt(C) replaced by C:

Compare these with the results using sqrt(C) :

Alan

[FanCG]

thanks a lot, Alan,

it is obvious that the latter one (1st order reaction scheme for k1) can much better fit the raw data with your codes.

while, the data I provide here is only one of 80 sets of data in my experiments.

in last year, we assumed that dA/dt=k1*C. but this postulate could not work better than the "half order reaction postulate" in a lot of other experimental conditions except for a few cases. so at last I adopted the half order assumption.

but this conclusion was drawn from the calculation in matlab and maple, I think I would compare two reaction schemes using remained raw data in mathcad according to your suggestion. If 1st order scheme can work better, it is also good news for myself, a lot of work can be simplified.

[jeanGiraud]

>my postulates on the reaction scheme are that ... <

______________________

It should be easy for Mathcad 11, Maybe not 15 ?

I don't understand sqrt(C) ? What's that doing

in a DE. It should write c = sqrt(C)

Typical reaction solved Mathcad 11

[FanCG]

the result is perfect while the reaction scheme is different from what i considered as following (for there won't be reverse reactions in my system, or at least no need to be considered )

according to around 80-90 sets of experiment data, the rough processes being carried out in the reactor are

gas C + solid_mixture == gas A

(half order reaction <k1> with respect to gas C)

dA/dt = k1 *m* sqrt(C)

gas C + solid_mixture == gas B

(first order reaction <k2> with respect to gas C)

dB/dt = k2*m*C

gas C (ads) + solid_mixture == C*||solid ==.....

(maybe first order adsorption <k3> with respect to gas C)

-dC(ads)/dt = k3 * m*C

hence consumption of C can be written as

-dC/dt = k1 *m* sqrt(C) + (k2+k3)*m*C

where m is a parameter concerned with solids bulk density, since i have got A(t), B(t) and C(t), so what to do next is using the data to fit the functions determined by the 3 blue equations.

[jeanGiraud]

Follow me Fan,

The magenta and blue fit have same model [their data sets make sense]. The rational fit [the one not done] does not make sense vs the magnitude of the coefficients but easy to fit too [ expect before midnight or sooner]. We can probably fit the three in one shot vis "MinerrOdesolve" ... but the MinerrOdesolve will not output the 3 model solutions, it will output only the data-set-solution. Still very interesting the data-set-solution whereas it will confirm the model DE system is correct. When a MinerrOdesolve does not fit at all, it indicates either the data collection is incorrect and/or the model DE system is incorrect. If you see my point, the MinerrOdesolve is like a cross correlation.

jmG

[jeanGiraud]

>See attached, ......... as M11 doesn't find a solution.<

_______________________

It does for me !

jmG

[FanCG]

it seems better than what i got in matlab.

do you use the same method as Alan's in ur mathcad 11?

[jeanGiraud]

Fan CG wrote:

it seems better than what i got in matlab.

do you use the same method as Alan's in ur mathcad 11?

We will soon get lost because there is confusion. The two fits are direct fits from the data sets that I manually composed. The third one wasn't included because the initialing parameters didn't work immediately, but that is not a conclusion. The 3 fits above are direct fit from the given DE [another Mathcad 11 project 2007].

jmG

[jeanGiraud]

Fan CG wrote:

it seems better than what i got in matlab.

do you use the same method as Alan's in ur mathcad 11?

Tomorrow, I will try the MinerrOdesolve". Again that will not give a function, but will enable checking if the DE system is correct. If the system is correctly dressed [modelled], it should fit the points as nicely as the proposed models here. You can reply in this thread, but the work sheet is in a more accessible conference "Usage" Fitting models [reaction rate]. You will find two models of extraordinary quality. The "generalized exponential" model is for you to evaluate if the first point of the decay is in doubt and could eventually be better estimated [ I have doubt on that point]. The second model "Transmute" fits the first point. Up to you again to evaluate if this model can fit again in the case of better measured value or initial estimate. You mention Alan model, but there is no work sheet ! Done in Mathcad 11, no sweat.

jmG

[FanCG]

thanks, jean

Alan's code is in this page at present, maybe you didnot notice it

In addition the 'function' i mentioned above is actually a kind of implicit relationship determined by the ODEs, not classic explicit functions.

BTW, the component C in the equations is the oxygen level in the reactor, the reproducibility of it can be considered good in different operating conditions.

following are some other simulation results i previously finished with matlab under the same postulates (the data is different from I posted in this thread, they have different temperature or solids mass or moisture)

it seems that the reaction scheme keeps working well (pls forgive me that I couldnt say too much details on the systems here, and can only post a small part of raw data).

condition 1

from top to bottom A,B,C

condition 2

from top to bottom : A,B,C

At last, I have to say that English is not my mother tongue, and there is a sentence in your reply i cannot get what you mean:

I want to know what you mean by "a more accessible conference "Usage". what make me confused is the word 'conference usage',

thank you for your tolerating my poor English

[AlanStevens]

jean,

Here is the worksheet I did saved in M11 format so you can see the method used. I haven't tested this particular worksheet in M11. ( Note that the raw data for C does not show on the first plot because 11 doesn't have a secondary y-axis, which is where Fan CG plotted it.)

I take your comments in the posts above to mean that you fitted some explicit function (such as a rational polynomial or something) to the raw data. If so, this doesn't address how to determine the rate constants for use in the ODE based model, which is what Fan CG wants.

Alan

[AlanStevens]

jean,

Ok. I've now looked at your reaction rate worksheet and you are indeed fitting the data to a couple of arbitrary explicit functions. My understanding (which might well be wrong) is that Fan CG is trying to develop a model based a little more on fundamental principles, using conservation laws, and is seeking to determine some of the constants of such a model by fitting the model to the data. I guess he then hopes to use the model to make predictions for situations for which he has no empirical data.

Alan

[ViktorKorobov]

Dear Fan CG!
In itself procedure of an estimation of rate constants with use Mathcad is not difficult. Here it is possible to apply various ways. Here I result a way with application of two solving blocks - Given/Odesolve and Given/Minerr. Construction of adequate kinetic model of process is much more more a challenge. Here I agree with Alan who has fairly noticed that the model with proportionality to a square root component C is inadequate to experiment. It is necessary to consider also an error of registration of experimental data. In the attached files I have analysed 4 cases. Your experimental data testify that model 1 is most preferable. This model assumes the first order of all reactions. Besides the model 1 considers an error of measurement of initial concentration of a component C.
Viktor

[FanCG]

thanks for your so many solutions, I will compare the 1st order postulate with the half order one using all of my remained data to check which can better describe the reaction

thanks

[jeanGiraud]

Alan,

You are correct and Viktor too about the DE underlying system. Most of the times we just have data from non expert originator and thus collaboration ends as a fit . An immense amount of work was done in the Category "Dose response", for Walter. I will put more of that stuff in this forum... just a matter of making a new collection. At least, this forum starts reviving with good material and fair visitors. Would you be kind enough to repost Viktor work sheet that belongs to your project, I think it should be possible to make it work in 11 from 15 [that will save me some precious time]. From what we have discussed with Walter, 15 takes the 11 setup.

Jean

[AlanStevens]

jean,

Attached are Viktor's worksheets in M11 format.

Alan

[ViktorKorobov]

Thanks Allan for conversion of files. Here I wish to make the small remark. I think always that MC11 is the best version. Here we have not numerous case of advantage of version 14. Solver AdamsBDF proves to be very effective in return inverse kinetics problem..

[jeanGiraud]

In other words of what you say that Mathcad 11 is the best version, sure : PTC Mathcad are only clones or cracks. Now the bad news is to come for you two wonderful collabs [Alan, Viktor] in this project. I have quickly examined Alan converted version, it opens and wants to run, it did run once and I could see the plots. But THEN appeared the malicious bug. What is the malicious bug ? For years, collabs had converted their 14, and more recently 15 to 11. What was convertible was fine and I could repost the corrected work sheet. But since few months, as I open the converted work sheets, it runs once and then either Mathcad 11 crashes or "an internal error has occurred" . That's a new bug in 14, 15 ... if it didn't exist before a certain date and started existing after, that's because of what I call "malicious bug". It seems to have appeared during the construction period of the "PTC Mathcad Facebook". Not surprising, this web community has been "spamed" [maybe still is], why not a virus that no AV will find. PTC may very well be forced to burn all their boxes and return to Mathcad 11.

All that to say that the poor MinerrOdesolve reveals the DE system is incorrect. It can't fit worse than a model fit. I will do my best, constructing the work sheet from blank. I have already done the data set. The fit from the "Transmute model" seems to indicate that in the DE some coefficients aren't simply numbers.

Cheers , Jean

[FanCG]

But why can't PTC do better in compatibilities among different versions of mathcad ?

it is said that the program and symbolic math would be cancelled in mathcad prime

[AlanStevens]

jean,

"Now the bad news is to come for you two wonderful collabs [Alan, Viktor] in this project. I have quickly examined Alan converted version, it opens and wants to run, it did run once and I could see the plots. But THEN appeared the malicious bug. What is the malicious bug ? For years, collabs had converted their 14, and more recently 15 to 11. What was convertible was fine and I could repost the corrected work sheet. But since few months, as I open the converted work sheets, it runs once and then either Mathcad 11 crashes or "an internal error has occurred" . That's a new bug in 14, 15 ..."

I've just downloaded the M11 converted format Viktor files and tested them in M11. I don't see the behaviour you report. They run ok in M11 (though they are just short of 15 times slower than in M14) and do not cause my computer or even Mathcad to crash. Did you just run them once, or did the crash happen after you modified something or what?

Alan

[jeanGiraud]

>I've just downloaded the M11 converted format Viktor files and tested them in M11. I don't see the behaviour you report. They run ok in M11 (though they are just short of 15 times slower than in M14) and do not cause my computer or even Mathcad to crash. Did you just run them once, or did the crash happen after you modified something or what ?<

PTC Mathcad versions are just not Mathcad [piece of scrap]. If it runs, that's one thing. If it crashes Mathcad, eventually the all box, or other kind of error , it results from the impossibility to to do anything significant in the converted work sheet. My observations are not casual and only pertaining to this particular transmission, no: my observations are anterior to a very recent from [ ...] that PTC had no other interest in Mathcad than their version 15 [which I call version sub-zero] ... This compobox does not even print on the screen as fast as I type 2 fingers, what would it be with 10 fingers !!!

All that to say : thanks for this wonderful project, all done in < 5 minutes "Glorious Mathcad 11". Calculation time < 5 seconds as long as the initials will initialize. How many seconds in your 14 ? Your DE system is as expected not correctly dressed.

Not sure what ERR means in that case ???

Jean

[FanCG]

yes,Alan, you got what I mean, the reaction rate constants are what I want at first