I have defined the function functiontoplot using the following piece of code:
f1 = c*(S*A^2)/H - mu*A + ρa*Y;
f2 = c*S*A^2 - ν*H + ρh*Y;
f3 = c0 - γ*S - ϵ*Y*S;
f4 = d*A - e*Y + Y^2/(1 + f*Y^2);
eq = Solve[{f1 == 0, f2 == 0, f3 == 0, f4 == 0}, {A, H, S, Y}][[3]];
jacobianmat = D[{f1, f2, f3, f4}, {{A, H, S, Y}}];
diffmatrix = DiagonalMatrix[{D1, D2, D3, D4}];
det = Simplify[Det[jacobianmat - μ*diffmatrix]];
der = Simplify[D[det, μ]];
ksquared = Solve[der == 0, μ][[1]];
c = 0.002;
mu = 0.16;
ρa = 0.005;
c0 = 0.02;
γ = 0.02;
ϵ = 0.1;
d = 0.008;
e = 0.1;
f = 10.0;
D1 = 0.001;
D2 = 0.02;
D3 = 0.01;
D4 = 10^-7;
det = Simplify[det];
ksquared = Simplify[ksquared];
eq = FullSimplify[eq];
functiontoplot = det /. ksquared /. eq;
functiontoplot depends on two parameters $\rho h$ and $\nu$. I know that it has a root approximately at the point where $\rho h=1.0011\cdot 10^{-5}$ and $\nu=0.04$. I've tried to plot the function using the following command:
ContourPlot[
functiontoplot == 0, {ρh, -0.001, 0.001}, {ν, 0.03, 0.05}]
However, after waiting for a long time, I get an empty plot as an output. I know that plotting an implicit function is an expensive process, so I want to implement an accurate faster method to plot my function.
Until now, I have tried with the pseudo-arclength continuation method explained in the first answer here:
In particular, after defining the function TrackRootPAL, I use the following commands:
tr = TrackRootPAL[{functiontoplot}, {ν}, {ρh, -0.001, 0.001},
1.0011*10^-5, {0.04}];
Plot[Evaluate[ν[ρh] /. tr], {ρh, -0.001, 0.001}]
Nevertheless, when it improves the initial condition, I get the complex number ν[0]==0.0399995 +9.22304*10^-21 I, so NDSolve is not able to keep working with the equation and, since it is a differential algebraic equation, there are no more methods available for NDSolve. I have even tried applying Chop to the initial condition but that then NDSolve needs to work with other complex numbers so I get warnings all of the time. Is there any solution for this using my approach or maybe another method to plot the curve I want to plot?
Subscript: https://mathematica.stackexchange.com/a/18395, point 3. Also: "Really the best way to use subscripts is to not use them." – Michael E2 Aug 12 '21 at 18:07functiontoplothas aLeafCountaround $5\times10^7$ and took over 90 sec. to evaluate once. -- And havingcandSubscript[c, 0]is a bad idea (see prev. comment). Probably not a big issue here, but I could be wrong. – Michael E2 Aug 12 '21 at 18:47functiontoplotthrough several functions, each computing a numerical step, instead of composing all of the symbolic expressions into a giant expression forfunctiontoplot. I mean the number of lines of code above is fairly short, the twoSolvebeing the only potentially long processes. But they're usually faster with fewer symbolic parameters. Something to consider if that's how the expressions get so large. – Michael E2 Aug 12 '21 at 18:56functiontoplotand now theLeafCountoffunctiontoplotis approximately equal to $8\times 10^6$. When you talked about refactoringfunctiontoplot, what do you mean with CS lingo? – edgardeitor Aug 12 '21 at 19:30eqhas 7 solutions andksquaredhas 3. How do you know that theeq[[3]]andksquared[[1]]are the ones you want? Are there any constraints on{A, H, S, Y, \[Mu]}? – Chris K Aug 12 '21 at 20:15eqrefers to the equilibrium points of the system. I only want to work with the stable one for the parameter values I am considering at the beginning. The equilibrium number three is the only one that fulfills that condition. On the other hand, for the case ofksquared, I am taking into account the solution that has lets $\rho_h=1.0011*10^{-5}$ and $\nu=0.04$ to be a common solution to the equations $det=0$ and $der=0$. I am following that approach because I need to know the value of $\mu$ at each of the points of the curve I need to plot. – edgardeitor Aug 12 '21 at 20:39{A, H, S, Y, \[Mu]}need to be real? positive? something else? – Chris K Aug 12 '21 at 21:38