0
$\begingroup$

I am solving the following ODE in Mathematica using NDSolve:

$$-dC_A/dt=k_sC_A^3+k_eC_A^2$$

where $k_s$ and $k_e$ are constants. The code I am using to solve the ODE is

s0 = NDSolve[{a'[t] == -ksc*a[t]^3 - kec*ccat*a[t]^2, a[0] == 2.5}, 
a, {t, 0, 2000}];

However, I have noticed that I get different results depending on when I end my time range. Specifically, I am solving

$$500=2.5/C_A(t)$$

for t.

If my time range goes to 2000, I get $t=4334.91$, but if my time range goes to 20,000, I get $t=17,705$. Any idea why this is occurring and which of the two is the correct result?

Edit: I've attached the code used to find t and parameter values.

NSolve[2.5/Evaluate[a[t] /. s0] == 500, t]

$k_s=0.2$ and $k_e=0.00825$

$\endgroup$
  • $\begingroup$ I'll assume a[t] in your code stands for $C_A$. How do you solve $500=2.5/C_A(t)$? If your maximum $t=2000$ in NDSolve, anything further would be extrapolation. $\endgroup$ – Chris K Apr 27 '17 at 18:49
  • $\begingroup$ @ChrisK You assumption is correct. Sorry if I was not clear, but in the first case it goes to 2000, but in the second case I've changed the code to $t_{\text{end}}=20,000$. $\endgroup$ – Tunk Apr 27 '17 at 18:51
  • $\begingroup$ Thanks. How do you solve $500=2.5/C_A(t)$? Do you work on the InterpolatingFunction after the NDSolve or use a WhenEvent? Have you tried plotting the solution to see which answer is correct? $\endgroup$ – Chris K Apr 27 '17 at 18:55
  • $\begingroup$ @ChrisK I've attached the code I'm using to solve for t. In regards to using a plot to determine which is the correct value, the plot also shows that the solution changes depending on the time scale. $\endgroup$ – Tunk Apr 27 '17 at 18:58
  • $\begingroup$ Hard to say for sure without knowing the parameter values, but I'm pretty sure that the real answer occurs at $t=17705$ and the other one is based on extrapolation of your InterpolatingFunction outside the range it was defined on. $\endgroup$ – Chris K Apr 27 '17 at 19:01
0
$\begingroup$

When you solve the differential equations to tmax=2000, the only way s0 can go beyond that time is by extrapolation, which easily goes wild. Look at this, where the valid solution is green and the extrapolated one is red:

ksc = 0.2;
kec = 0.00825;
ccat = 1;
tmax = 2000;
s0 = NDSolve[{a'[t] == -ksc*a[t]^3 - kec*ccat*a[t]^2, a[0] == 2.5}, a, {t, 0, tmax}];
Show[
  Plot[Evaluate[{500, 2.5/a[t] /. s0}], {t, 0, 2000}, PlotStyle -> {Blue, Green}],
  Plot[Evaluate[{500, 2.5/a[t] /. s0}], {t, 2000, 20000}, PlotStyle -> {Blue, Red}],
  PlotRange -> {{0, 20000}, {0, 600}}
 ] 

Mathematica graphics

Here we actually solve the differential equations for tmax=20000, which means no extrapolation is used:

tmax = 20000;
s0 = NDSolve[{a'[t] == -ksc*a[t]^3 - kec*ccat*a[t]^2, a[0] == 2.5}, a, {t, 0, tmax}];
Plot[Evaluate[{500, 2.5/a[t] /. s0}], {t, 0, tmax}, PlotRange -> All, PlotStyle -> {Blue, Green}]

Mathematica graphics

Clearly the real answer is around $t=17700$.

| improve this answer | |
$\endgroup$
  • $\begingroup$ Thanks a ton! Completely makes sense now. $\endgroup$ – Tunk Apr 27 '17 at 19:22

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.