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:

 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:

Is there any predictor-corrector method in Mathematica for solving nonlinear system of algebraic equations?

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?

  • $\begingroup$ FYI, re Subscript: mathematica.stackexchange.com/a/18395, point 3. Also: "Really the best way to use subscripts is to not use them." $\endgroup$
    – Michael E2
    Aug 12 at 18:07
  • $\begingroup$ functiontoplot has a LeafCount around $5\times10^7$ and took over 90 sec. to evaluate once. -- And having c and Subscript[c, 0] is a bad idea (see prev. comment). Probably not a big issue here, but I could be wrong. $\endgroup$
    – Michael E2
    Aug 12 at 18:47
  • $\begingroup$ Yeah maybe you are right! I removed all the subscripts from the script. However, the result is still the same. $\endgroup$ Aug 12 at 18:48
  • $\begingroup$ Sorry, but I don't have anything really helpful to suggest. Possibly you could refactor (in CS lingo) functiontoplot through several functions, each computing a numerical step, instead of composing all of the symbolic expressions into a giant expression for functiontoplot. I mean the number of lines of code above is fairly short, the two Solve being 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. $\endgroup$
    – Michael E2
    Aug 12 at 18:56
  • 1
    $\begingroup$ For your reference parameters, eq has 7 solutions and ksquared has 3. How do you know that the eq[[3]] and ksquared[[1]] are the ones you want? Are there any constraints on {A, H, S, Y, \[Mu]}? $\endgroup$
    – Chris K
    Aug 12 at 20:15

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.