Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

Good morning,

I have an ODE linear system $\dot{\vec f} = A \vec f$ depending parametrically to two real parameters, s and p. I'd like to get a density map of the solutions of $f_A/(f_A+f_B)$ at $t=20$, in the space of parameters ${s, p} \in [0; 1]$.

I post my try:

rA = 2.1
rB = 1.0
fA := ParametricNDSolveValue[{A'[ t] == (rA - 1 + p) A[ t] + s B[t], B'[t] == (1 - p) A[t] + (rB - s) B[t], A[0] == 2, B[0] == 2}, A, {t, 0, 20}, {{s, p} \[Element] Reals}]
fB := ParametricNDSolveValue[{A'[ t] == (rA - 1 + p) A[ t] + s B[t], B'[t] == (1 - p) A[t] + (rB - s) B[t], A[0] == 2, B[0] == 2}, B, {t, 0, 20}, {{s, p} \[Element] Reals}]
DensityPlot[fA[s][p]/(fA[s][p] + fB[s][p]), {s, 0, 1}, {p, 0, 1}]

Here there's what I get

The plot I obtain can't be what I need, as I never ask Mathematica to evaluate the functions at $t=20$. I'm really new to Mathematica and I've the impression to be missing something really basic about its semantics. I hope to learn from your help, thank you in advance!

share|improve this question
    
I know that nowhere in my code there is $t=20$, but when I tried few things to do that I got errors like "ParametricNDSolve::dsvar: 20 cannot be used as a variable" –  Orso Jul 8 at 9:06
1  
I cut and pasted your code and the plot works? Try Clear[t] first, it looks like you may have executed t=20 before. –  Ymareth Jul 8 at 9:16
    
Hi @Ymareth and thank you for your help. Indeed, I already do Clear[t] before the lines I posted but it doesn't change anything in my case... –  Orso Jul 8 at 9:21
    
You're right, I pasted it into a new notebook and I obtain a plot. The problem is that nobody told the machine to plot at $t=20$, so I guess the plot is not at all what I need... –  Orso Jul 8 at 9:38

1 Answer 1

up vote 1 down vote accepted

I may have misunderstood the aim. I post this as motivating:

I have re-written the system of equations (no guarantee error free):

mat = {{1.1 + p, s}, {1 - p, 1 - s}}; pf = 
 ParametricNDSolveValue[{{a'[t], b'[t]} == 
    mat.{a[t], b[t]}, {a[0], b[0]} == {2, 2}}, {a[t], b[t]}, {t, 0, 
   20}, {s, p}];
foi[par1_, par2_, u_] := (#1/(#2 + #1) & @@ pf[par1, par2]) /. t -> u

Looking at the solutions $a(t),b(t)$:

Manipulate[
 Column[{ParametricPlot[pf[v1, v2] /. t -> w, {w, 0, 20}, 
    PlotRange -> {{0, 10^12}, {0, 10^12}}, Frame -> True, 
    FrameLabel -> {"a", "b"}], 
   LogPlot[Evaluate[pf[v1, v2] /. t -> w], {w, 0, 20}, 
    PlotLegends -> {"a(t)", "b(t)"}
    , PlotRange -> {0, 10^12}, Frame -> True, 
    FrameLabel -> {"t", "a|b"}]}], {{v1, 0, "s"}, 0, 
  1}, {{v2, 0, "p"}, 0, 1}]

enter image description here

As expected with these linear system with constant positive coefficients->exponentials

The desired density plot for function (foi evaluated at t=20) above:

DensityPlot[foi[x, y, 20], {x, 0, 1}, {y, 0, 1}, 
 ColorFunction -> "SunsetColors", PlotLegends -> Automatic, 
 MeshFunctions -> (#3 &), Mesh -> {Range[0, 1, 0.1]}, 
 FrameLabel -> {"s", "p"}]

enter image description here

or in 3D:

Plot3D[foi[x, y, 20], {x, 0, 1}, {y, 0, 1}, 
 ColorFunction -> (ColorData["SunsetColors"][#3] &), 
 MeshFunctions -> (#3 &), Mesh -> {Range[0, 1, 0.1]}, 
 AxesLabel -> {"s", "p", 
   "\!\(\*FractionBox[\(a \((20)\)\), \(a \((20)\) + b \
\((20)\)\)]\)"}]

enter image description here

Apologies if I have misinterpreted your aim.

share|improve this answer
    
Exactly, @ubpqdn, that's what I wanted! Thank you so much, I've learnt a lot from your solution! –  Orso Jul 8 at 13:48

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

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