2
$\begingroup$

I have a system of coupled PDEs with five parameters that I only know ranges for rather than explicit values, and I was wondering if I could define a parametric equation with ParametricNDSolve and then use FindFit with some data to estimate these parameters? My problem is that I have some coupled equations.

Firstly I have a set up equation which seems to work fine, given by 

(*Constants!*)
a = 10*10^-7; omega = 3.0138*10^7; Do2 = 2*10^-9; po = 106;
ro = 263*10^-6; micron = 
 1*10^-6; k = 1; (*eo = 100;*)(* qm = 10^-4;*)
mic = 1*10^-6; De = 5.5*10^-11;

(*First equation - set up eq!*)
s = Quiet[
   NDSolve[{D[Ox[r, t], t] - 
       Do2*(D[Ox[r, t], r, r] + (2/r)*(D[Ox[r, t], r])) + (a*
          omega)*((Ox[r, t])/(Ox[r, t] + k)) == 0, Ox[r, 0] == 0, 
     Ox[micron, t] == 0, Ox[ro, t] == po}, 
    Ox, {r, micron, ro}, {t, 0, 14400}]];
p = Ox /. First[s];

Which gives me a function p[r,t] . I then use this in a parametric set up with parameters $eo,qm,kme,kmn$ and $j$ and initial / BC conditions to solve parametrically for $Ef1$;

(*Second equation*)
eqnDe = D[Ef1[r, t], t] - 
   De*(D[Ef1[r, t], r, r] + (2/r)*(D[Ef1[r, t], r])) + 
   qm*((kme)/(kme + 
      p[r, t])*((p[r, t])/(p[r, t] + 
       kmn)) + (1 - (p[r, t])/(p[r, t] + kmn))*j)*Ef1[r, t];

l = ParametricNDSolve[{eqnDe == 0, Ef1[r, 0] == 0, 
Derivative[1, 0][Ef1][micron, t] == 0, Ef1[ro, t] == eo}, 
   Ef1, {r, micron, ro}, {t, 0, 14400}, {kme, kmn, j, qm, eo}];

b = Ef1 /. l;

This seems to work so far (though Im unfamiliar with ParametricNDSolve) but my problem is that I want to use the output of this and the parameters to find a new value Eb1, related to the above solution by;

$\frac{dEb1}{dt} = qm\left(\frac{p}{p + kmn}\frac{kme}{p + kme} + (1 - \frac{p}{p + kmn})j \right)Ef1$

with the initial condition $Eb1[r, 0] == 0$ - the problem is I have tried coded this in and have had no luck. I tried

eqnBo = D[Eb1[r, t], 
t] - (b[r, 
  t])*(((p[r, t])/(p[r, t] + kmn))*((kme)/(kme + 
       p[r, t])) + (1 - (p[r, t])/(kmn + p[r, t]))*j);

 x = ParametricNDSolve[{eqnBo == 0, Eb1[r, 0] == 0}, 
    Eb1, {r, micron, ro}, {t, 0, 14400}, {kme, kmn, j, qm, eo}];

And this produces an error;

ParametricNDSolve::fpct: "Too many parameters in {kme,kmn,j,qm,eo} to be filled from {r,t}"

Anyone have any idea if what I want to do is possible, and whether this error simply means there's too much to fit? The data I want to fit it to is here.

I'd be really grateful for any guidance I can get on this, or improved methods of finding the data I want!

Improve this question
$\endgroup$
5
  • $\begingroup$ I should also point out similar error message EVEN if you specify every parameter bar one; for example if I specify j, kmn, eo and qm, with only parameter kme I get; ParametricNDSolve::fpct: Too many parameters in {kme} to be filled from {r,t} $\endgroup$
    – DRG
    Commented Mar 7, 2014 at 15:58
  • 1
    $\begingroup$ There seems a discordance between Ox[r, 0] == 0 and Ox[ro, t] == po in your boundary conditions. $\endgroup$
    – qwerty
    Commented Mar 8, 2014 at 3:01
  • 1
    $\begingroup$ One problem I see is that you are using the your parametric function b incorrectly; it takes on the form b[kme, kmn, j, qm, eo] and you are using b[r,t]. $\endgroup$
    – bobthechemist
    Commented Mar 8, 2014 at 13:56
  • $\begingroup$ Two clarifications: (1) It looks like p reaches a steady state around t = 10. Is there a reason why t -> 14400 is needed in the analysis? (2) Can you take another look at the data you posted? It is only one column and it's not clear to me what it is. $\endgroup$
    – bobthechemist
    Commented Mar 8, 2014 at 14:47
  • $\begingroup$ Thanks guys - You're absolutely right re: b function - stupid oversight on my part! p is essentially steady state, and even though I have a analytical function for this case, so I may implement that. Regarding espirit's observation - this is absolutely correct but I'm not sure how to work around it. I've expanded on this under Michael E2's excellent answer. Thanks again folks - really helpful for a newb! $\endgroup$
    – DRG
    Commented Mar 9, 2014 at 18:10

1 Answer 1

Reset to default
3
$\begingroup$

If you avoid the intermediate integration of Ef1, you can do it all at once:

eqnBo = D[Eb1[r, t], t] - (Ef1[r, t]) * (((p[r, t])/(p[r, t] + kmn)) * 
          ((kme)/(kme + p[r, t])) + (1 - (p[r, t])/(kmn + p[r, t]))*j);

x = ParametricNDSolve[{eqnBo == 0, Eb1[r, 0] == 0, eqnDe == 0, 
    Ef1[r, 0] == 0, Derivative[1, 0][Ef1][micron, t] == 0, 
    Ef1[ro, t] == eo}, 
   Eb1, {r, micron, ro}, {t, 0, 14400}, {kme, kmn, j, qm, eo}];

I still get the error message @esprit alludes to.

Eb1[0.1, 0.1, 0.1, 0.1, 0.1] /. x

ParametricNDSolve::ibcinc: Warning: boundary and initial conditions are inconsistent. >>

(* InterpolatingFunction[{{1. 10^-6, 0.000263}, {0., 14400.}}, <>] *)
Improve this answer
$\endgroup$
1
  • $\begingroup$ Excellent answer, cheers - I have kind of ignored this as I've found it difficult to code up; in the region $o \leq r \leq r_{o}$ , $p[r,0] = 0$. However, beyond $r_{o}$, $p[r_{beyond},t] = p_{o}$ - I have no idea how to code this interface, so I think this leads to the issue of conflict. If anyone knows how I tell Mathematica that $p_{o}$ is the value BEYOND $r_{o}$ I'd be extremely grateful to learn syntax! $\endgroup$
    – DRG
    Commented Mar 9, 2014 at 18:13

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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