2
$\begingroup$

I have the following ODE that I have solved with NDSolve:

y'[x] == -l*a*(y[x]^2 - c*y[x] - p*x^3*Exp[-2*x])/x^2

where $l = 1.508*10^{13}$, $a = 5*10^{-9}$, $c = 10^{-11}$, and $p = 1.038*10^{-5}$ are constants using the code 

sol1 = NDSolve[{y'[x] == -l*a*(y[x]^2 - c*y[x] - p*x^3*Exp[-2*x])/x^2,
y[10] == 10^(-5)}, y, {x, 10, 100}, MaxStepSize -> 0.001].

But now that I've done this, I want to vary $a$. Specifically, I want to graph $y[100]$ as a function of $a$, where $a$ varies from $10^{-7}$ to $10^{-11}$.

I've tried looking at past threads but I haven't seen anyone do what I want to do.

I've tried making sol1 a function of $a$ but that doesn't work because then I've fed NDSolve a non-numerical value.

I've tried solving the ODE as a PDE (not an ideal solution...) by taking $a$ as a variable with the following code

sol2 = 
NDSolve[{D[y[x, a], x] == -l*a*(y[x, a]^2 - c*y[x, a] - p*x^3*Exp[-2*x])/x^2, 
y[10, a] == 10^(-5)}, y, {x, 10, 100}, {a, 10^-11, 10^-7}, MaxStepSize -> 0.001].

However, I run into an error anyway: "At x == 10.`, step size is effectively zero; singularity or stiff system suspected."

Does anyone know how to graph the end result of sol1, $y[100]$, as a function of $a$?

Improve this question
$\endgroup$
7
  • $\begingroup$ If NDSolve runs fast in your case, you can solve the equation multiple times for differene $a$, retrieve the value of $y(100)$ for each and plot a ListPlot. $\endgroup$
    – corey979
    Commented Sep 7, 2016 at 18:15
  • $\begingroup$ I could, but it'd really be ideal to get a continuous function. I'm trying to replicate the results of a paper, and their $f[100]$ v $a$ plot is continuous. $\endgroup$
    – Ecclesiastic
    Commented Sep 7, 2016 at 18:21
  • $\begingroup$ But you know that in a computer nothing is contiuous? Even when do Plot[x^2,{x,0,1}] you in fact obtain a set of connected points, not a mathematically continuous graph. You can plot yours for e.g. 100 or 1000 values of $a$. $\endgroup$
    – corey979
    Commented Sep 7, 2016 at 18:25
  • $\begingroup$ I suppose. But I think the difference would be that if I want to obtain $y[100]$ for, say, $a = 1.5*10^{-11}$ and I only have calculated $y[100]$ for the points $a = 10^{-11}, 2*10^{-11}, ..., 10^{-7}$ beforehand, then I will have to go calculate $y[100]$ again, whereas if I want the value of $x^2$ for any $x$, I can obtain it much more easily. So that's the sort of sense in which I'd like the function to be continuous. It seems like that should be possible. $\endgroup$
    – Ecclesiastic
    Commented Sep 7, 2016 at 18:33
  • 3
    $\begingroup$ ParametricNDSolve often is used for this purpose. $\endgroup$
    – bbgodfrey
    Commented Sep 7, 2016 at 19:07

1 Answer 1

Reset to default
3
$\begingroup$

As mentioned by bbgodfrey in his comment, you can use ParametricNDSolve, or otherwise ParametricNDSolveValue:

l = 1.508 10^13;
c = 10^(-11);
p = 1.038 10^(-5);

eqs = {y'[x] == -l*a*(y[x]^2 - c*y[x] - p*x^3*Exp[-2*x])/x^2, y[10] == 10^(-5)};

sol = ParametricNDSolveValue[eqs, y, {x, 10, 100}, {a}]

enter image description here

Plot[sol[a][100], {a, 10^-11, 10^-7}]

enter image description here

Improve this answer
$\endgroup$
1
  • $\begingroup$ Thanks so much! This is what I needed. As a quick note, I actually copied the l value wrong--it's $3.01593*10^{21}$. This results in bad behavior for the above graph but this can be fixed by adding AccuracyGoal -> 20 in the ParametricNDSolveValue. $\endgroup$
    – Ecclesiastic
    Commented Sep 8, 2016 at 5:19

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.