0
$\begingroup$

Edit: My approach was wrong, I have solved this - please see comments

I've been scouring the answers and the docs but I'm having trouble figuring out the answer to this one. In particular, Olksander R.'s answer was particularly helpful.

I have a single PDE I'm trying to fit to a dataset (R is a constant):

 a'[T] == A/B Exp[-Ea/(R T)] (1 - a[T])^n a[T]^m

Five total data points were collected, and they define more than one of the variables in the partial differential equation. Please note that the dadT variable represents the left-hand-side of the PDE.

BVal = {2, 5, 10, 15, 20};
Tval = {448, 461, 473, 480, 484};
aVals = {0.659090909, 0.617021277, 0.58, 0.568627451, 0.528301887};
dadT = {0.025, 0.023404255, 0.02, 0.021960784, 0.016981132};

Plugging these values in to the PDE gives me a set of 5 equations, with some of the variables parameterized and some needing to be fit.

I'm not sure how to handle the fact that I have data for the "left-hand-side" of the equation.

Edit: I added an equality of the a'[T] value with it's measured value.

I still intend to fit (Transpose[{aVals,dadT}]) with NonlinearModelFit later.

a'[T]==(A E^(-((0.120272 Ea)/T)) (1-a[T])^n a[T]^m)/B
a'[448]==0.025==1/2 0.340909^n 0.659091^m A E^(-0.000268465 Ea)
a'[461]==0.0234043==1/5 0.382979^n 0.617021^m A E^(-0.000260894 Ea)
a'[473]==0.02==1/10 0.42^n 0.58^m A E^(-0.000254275 Ea)
a'[480]==0.0219608==1/15 0.431373^n 0.568627^m A E^(-0.000250567 Ea)
a'[484]==0.0169811==1/20 0.471698^n 0.528302^m A E^(-0.000248496 Ea)

Where now I need to determine {A, B, n, m, Ea}.

I found some great answers using ParametricNDSolveValue to create parameterized system for NonlinearModelFit to then regress to the values of variables. Edit: My system of equations spits out a ParametricFunction, but it is not readily yielding and interpolating function to then fit:

nDss = ParametricNDSolveValue[system, a, {T, 448, 484}, {A, n, m, B, Ea}]
ParametricFunction[Expression: a, Parameters: {A,n,m,B,Ea}]

Attempting to parameterize the solution to get a single interpolating function gives me the following error:

nDss[.3, 2,.7, .8, .001]
ParametricNDSolveValue::icordinit: The initial values for all the dependent 
variables are not explicitly specified. NDSolve will attempt to find 
consistent initial conditions for all the variables.
ParametricNDSolveValue::ndnco: The number of constraints (5) (initial 
conditions) is not equal to the total differential order of the system plus 
the number of discrete variables (1).

If I understand, Mathematica is trying to tell me the I don't have enough inital conditions to solve the PDE. Is that right?

Really, I don't need a solution to the PDE; I need to do a multivariate least-squares fit of my dataset to determined the variables in the PDE. My thinking in using this approach was that my data collection would allow me to parameterize the PDE enough to get some interpolants which could then be used to minimize, like in Oleksandr's answer linked above. Maybe my approach is wrong - any comments are helpful!!!

$\endgroup$
  • $\begingroup$ Repeat what you posted using a[T] instead of a. What do you get? Also, to what parameter dadT corresponds to? Is it a'[T] ? $\endgroup$ – Anton Antonov Apr 6 '18 at 15:35
  • $\begingroup$ I get the same error with a[t]': ParametricNDSolveValue::dvnoarg: The function a appears with no arguments. dadT corresponds to the value a'[t] - the derivative of a vs. T $\endgroup$ – Jake M Apr 6 '18 at 16:31
  • 1
    $\begingroup$ Write a'[T] not a[T]' $\endgroup$ – m_goldberg Apr 7 '18 at 1:47
  • $\begingroup$ Sorry, I made a typo in my comment - I didn't mean to type a[T]', I meant to type a[T]. All three function variables a, a[T], and a'[T] give me the exact same ` function appears with no arguments` error. $\endgroup$ – Jake M Apr 8 '18 at 15:31
  • $\begingroup$ Replace a with a[T] in your equation a'[T] == A/B Exp[-Ea/(R T)] (1 - a[T])^n a[T]^m $\endgroup$ – swish Apr 10 '18 at 18:28
0
$\begingroup$

It seems I was wrong-headed in my approach. Mathematica was trying to hint that the system was overdetermined - there were more boundary conditions than the order of the equation, which leads to infinite solutions. This should have led me to understand that despite my expression being a PDE, no derivatives appear on the right-hand side of equation, and I can measure the derivative left-hand side directly, so I can regress the right-hand side as an expression directly to the dataset.

Simply using NonlinearModelFit or FindFit on the properly-formatted matrix of measured parameters gives the right answer.

The variables go first in the matrix, followed by the estimate of the expression, regardless if they have a measured datapoint or not, which is what I was confounded by.

data=Transpose[{aVals, BVal, Tval, dadT}];
data //TableForm
0.659091    2   448 0.025
0.617021    5   461 0.0234043
0.58    10  473 0.02
0.568627    15  480 0.0219608
0.528302    20  484 0.0169811

Now, it's as simple as making sure your variable are listed in the same order...

nlm = NonlinearModelFit[data,
 {A/B Exp[(-Ea/(8.314472 T))] (1 - a)^n a^m,
 {n > 0, A > 0, Ea > 0, m > 0}}, 
 {A, n, m, Ea},  {a, B, T}];

nlm["BestFitParameters"]
{A -> 312.008, n -> 6.55572, m -> 0.831376, Ea -> 4749.2}

The errors are horrendous, due to the fact that the model is assumed to be unconstrained and the dataset is somewhat weak. I've yet to figure out if the fitted result is meaningful; that's something for later.

nlm["RSquared"]
0.97805

nlm["ParameterTable"]
FittedModel::constr: The property values {ParameterTable} assume an 
unconstrained model. The results for these properties may not be valid, 
particularly if the fitted parameters are near a constraint boundary.

    Estimate    Standard Error  t-Statistic P-Value
A   312.008 9715.21 0.0321154   0.979562
n   6.55572 28.0232 0.233939    0.853701
m   0.831376    28.5168 0.0291539   0.981445
Ea  4749.2  185893. 0.025548    0.983739
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.