Plot the solution of a differential equation

Consider: $$\frac{1}{(1-e)^2 w}=v-z$$ where $$v=f(e,i,L,W,a,n)$$, $$z=g(e,i,L,W,a,w,n)$$, and $$L = h(e,w,a,b)$$ are given functions.

From the above equation, $$e$$ can be derived as a function of $$w$$ along with the other parameters, i.e., $$e=e(w;\dots)$$. Using the latter, we can consider the following differential equation: $$\frac{\partial e}{\partial w}=\frac{e}{w}$$ I would like find $$e$$ and $$w$$ that satisfies the above differential equation. Obviously, both $$e$$ and $$w$$ will come out as a function of $$a$$ along with the other parameters.

Finally, I would like to create four separate plots for $$e$$, $$w$$, $$\frac{e}{w}$$, and $$\frac{ae}{w}$$ each of them against $$a\in [0,1]$$ with the parameter values of $$n=1$$, $$W=60$$, $$i=0.1$$, and $$b=0.6$$.

By closely referring to Michael E2's answer to this post, I tried to come up with the following code:

Clear["Global*"];
n = 1; W = 60;
L = (w/(b (a e)^b))^(1/(b - 1));
v = -(((-1 + e - i) (-i + (-1 + e + e i) L n) (1/((-1 + e) (-1 + e - i) w) + w/(-1 + e - i) + ((-1 + e) (1 - e L n) (-(1 + i)^2 (-1 + (1 + i)^((1 - e L n)/(L n - e L n)))^2 + a^2 i^2 (-1 + (1 + i)^(1 + 1/(L n - e L n)))^2 W^2))/(a i (1 + i) (1 - e + i) (1 - (1 + i)^((1 - e L n)/(L n - e L n))) (1 - (1 + i)^(1 + 1/(L n - e L n))) (i - (-1 + e + e i) L n) W)))/(i (1 + i + L n + e (-1 + (-2 + e - i) L n))));
z = -((a i (1 + i)^((1 + e)/(1 - e)) (-1 + (1 + i)^((1 - e L n)/(L n - e L n))) (-1 + (1 + i)^(1 + 1/(L n - e L n))) L n W + a (-1 + e) i (1 + i)^((1 + e)/(1 - e)) (-1 + (1 + i)^((1 - e L n)/(L n - e L n))) (-1 + (1 + i)^(1 + 1/(L n - e L n))) L n w^2 W + (-1 + e - i) (1 + i)^(-((2 e)/(-1 + e))) (-1 + e L n) w (1 + 2 i + i^2 - a^2 i^2 W^2 - 2 (1 + i)^(1 + 1/(L n - e L n)) ((1 + i)^(2 + 1/(-1 + e)) - a^2 i^2 W^2) + (1 + i)^(2 + 2/(L n - e L n)) ((1 + i)^((2 e)/(-1 + e)) -a^2 i^2 W^2)))/(a i^2 (1 + i) ((1 + i)^(e/(1 - e)) - (1 + i)^(1/(L n - e L n))) ((1 + i)^(e/(1 - e)) - (1 + i)^((1 + L n)/(L n - e L n))) (1 + i + L n + e (-1 + (-2 + e - i) L n)) w W));

myEQ = 1/((1 - e)^2 w) == v - z;

Block[{n = 1, W = 60, i = 1/10, b = 6/10, a = 1/2}, {#, Dt[#, w]} &@myEQ] /. {e -> e[w]} /. {e'[w] -> e[w]/w} /. {e[w] -> e} // Simplify;
icEQ = % /. Equal -> Subtract // Simplify;
icSOL = NSolve[icEQ == 0 && 0 < e < 1 && w > 0, {e, w}]

{e, w} = NDSolveValue[{ode = Solve[Block[{n = 1, W = 60, i = 1/10, b = 6/10}, D[{myEQ /. Equal -> Subtract, {e, w} . D[myEQ /. Equal -> Subtract, {{e, w}}]} /. {w -> w[a], e -> e[a]}, a] == 0], {e'[a], w'[a]}] /. Rule -> Equal, e[1/2] == (e /. First@icSOL), w[1/2] == (w /. First@icSOL)}, {e, w}, {a, $MachineEpsilon, 1}] ListLinePlot[e, PlotRange -> All, AxesLabel -> {a, "e"}] ListLinePlot[w, PlotRange -> All, AxesLabel -> {a, "w"}]  The code runs forever. Also, I failed to come up with a code that creates the four separate plots. By the way, the code is exactly the same as Michael E2's except that the equation in MyEQ is now a more complicated one with v, z, and L. Since Michael E2's code with the original MyEQ equation works well, I don't think there is an error in the code itself. Any help will be greatly appreciated. • You assign a value to b but you don't use that variable. Is b actually beta? Commented May 23 at 2:04 • @BobHanlon, yes, you are correct. I have edited accordingly and now they are all in b consistently. Thanks! – ppp Commented May 23 at 13:31 • You still have beta in your code. But even if I replace it with b your Block construct is hanging apparently because of the Simplify. Removing the Simplify it quickly completes. Next, removing Simplify on icEQ runs quickly but then ICSOL starts hanging. – josh Commented May 23 at 16:26 • Another suggestion: IcEQ as written above is in brackets as in {icEQ}. Not sure that's a construct error. However, if I remove the bracket and pick a value of e and then use FindRoot to solve w at 25 digits of precision I can get a solution as in: FindRoot[( icEQ[[1]] == 0 /. e -> 1/5), {w, 1/2}, WorkingPrecision -> 25]  returning {e,w)={1/5, 3.504}. If this is usable, you could generate a table of such FindRoots for your needs. – josh Commented May 23 at 16:51 • @ppp It is not clear from your code what is the domain of parameter a ? Commented May 27 at 4:59 1 Answer We can solve this problem as an optimization problem as follows. First note, that we don't know upfront are there real solutions for a given parameters or not. Therefore we need to compute a global minimum for myEQ to define a real solution. Differential equation $$\frac{\partial e}{\partial w}=\frac{e}{w}$$ can be integrated in general form as $$w=c e$$, where c is arbitrary constant. The global minimum is defined in a form L = (w/(b (a e)^b))^(1/(b - 1)); v = -(((-1 + e - i) (-i + (-1 + e + e i) L n) (1/((-1 + e) (-1 + e - i) w) + w/(-1 + e - i) + ((-1 + e) (1 - e L n) (-(1 + i)^2 (-1 + (1 + i)^((1 - e L n)/(L n - e L n)))^2 + a^2 i^2 (-1 + (1 + i)^(1 + 1/(L n - e L n)))^2 W^2))/(a i (1 + i) (1 - e + i) (1 - (1 + i)^((1 - e L n)/(L n - e L n))) (1 - (1 + i)^(1 + 1/(L n - e L n))) (i - (-1 + e + e i) L n) W)))/(i (1 + i + L n + e (-1 + (-2 + e - i) L n)))); z = -((a i (1 + i)^((1 + e)/(1 - e)) (-1 + (1 + i)^((1 - e L n)/(L n - e L n))) (-1 + (1 + i)^(1 + 1/(L n - e L n))) L n W + a (-1 + e) i (1 + i)^((1 + e)/(1 - e)) (-1 + (1 + i)^((1 - e L n)/(L n - e L n))) (-1 + (1 + i)^(1 + 1/(L n - e L n))) L n w^2 W + (-1 + e - i) (1 + i)^(-((2 e)/(-1 + e))) (-1 + e L n) w (1 + 2 i + i^2 - a^2 i^2 W^2 - 2 (1 + i)^(1 + 1/(L n - e L n)) ((1 + i)^(2 + 1/(-1 + e)) - a^2 i^2 W^2) + (1 + i)^(2 + 2/(L n - e L n)) ((1 + i)^((2 e)/(-1 + e)) - a^2 i^2 W^2)))/(a i^2 (1 + i) ((1 + i)^(e/(1 - e)) - (1 + i)^(1/(L n - e L n))) ((1 + i)^(e/(1 - e)) - (1 + i)^((1 + L n)/(L n - e L n))) (1 + i + L n + e (-1 + (-2 + e - i) L n)) w W)); eq1 = 1/((1 - e)^2 w) - v + z; eq2 = w - c e; var = {a, e, b, i, n, W, c}; con = {0 < a < 1, 0 < e < 1, 0 < c e < 2, W > 0, i > 0, n > 0, b > 0}; sol = NMinimize[{eq1^2 /. w -> c e, con}, var] (*Out[]= {3.05276*10^-21, {a -> 0.12749, e -> 0.131433, b -> 0.154335, i -> 0.170289, n -> 0.183486, W -> 0.19822, c -> 0.193286}}*)  Using this data we can compute also a local minimum, we have s2 = Drop[sol[[2]], 2] (*Out[]= {b -> 0.154335, i -> 0.170289, n -> 0.183486, W -> 0.19822, c -> 0.193286}*) var1 = {a, e}; con1 = {0 < a < 1, 0 < e < 1}; ini1 = Table[{var1[[j]], var1[[j]] /. sol[[2]]}, {j, Length[var1]}]; sol2 = FindMinimum[({eq1^2, con1} /. w -> c e) /. s2, ini1] Out[]= {8.26643*10^-18, {a -> 0.417144, e -> 0.497054}}  The way from the global minimum to the local minimum is some curve $$e=e(a)$$, that we can extract from FindMinimum using StepMonitor pts = Reap[ FindMinimum[{eq1^2 /. w -> c e /. s2, con1}, ini1, StepMonitor :> Sow[{a, e, eq1 /. w -> c e /. s2}], MaxIterations -> 1000]][[2, 1]];  Finally we can plot pts with contours of expression eq1==0, we have plot1=ContourPlot[(eq1 /. w -> c e /. s2) == 0, {a, .001, .5}, {e, .01, .99}, Epilog -> {Red, Line[pts[[All, 1 ;; 2]]], Point[pts[[All, 1 ;; 2]]]}, PlotLegends -> Automatic, FrameLabel -> {"a", "e"}]  To plot all other curves we can extract list of points {a,e} from plot1, then construct lists for w and plot it lst = plot1[[1]][[1]][[1]]; = Transpose[{lst[[All, 1]], c lst[[All, 2]]}] /. sol[[2]]; ListLinePlot[lstw, Frame -> True, FrameLabel -> {"a", "w"}]  All other functions like w/e=c and a e/w=a/c can be plot as well for this kind of solutions. Update 1 We can use ContourPlot to visualize solution for a given parameters as follows n = 1; W = 60; i = 0.1; b = .6;Table[ContourPlot[(eq1 /. w -> c e) == 0, {a, 0.01, .99}, {e, .4, .99}, PlotLabel -> Row[{"c = ", c}], FrameLabel -> {"a", "e"}, PlotPoints -> 100], {c, .02, .22, .05}]  The bottom branch in the picture above can be reproduced with FindRoot as sol10 = Table[ Select[Quiet@ Table[{a, e, Abs[eq1 /. w -> c e]} /. Quiet@FindRoot[(eq1 /. w -> c e) == 0, {e, .9}, Method -> {"Newton", "StepControl" -> "TrustRegion"}], {a, .001, .999, .001}], #[[2]] > 0 && #[[2]] < 1 && #[[3]] < 10^-10 &], {c, .02, .22, .05}]; p1 = ListPlot[Table[sol10[[j]][[All, 1 ;; 2]], {j, 5}], PlotLegends -> Table[Row[{"c =", c}], {c, .02, .22, .05}], AxesLabel -> {"a", "e"}, PlotRange -> {{0, 1}, {.4, 1}}]; sol04 = Table[ Select[Quiet@ Table[{a, e, Abs[eq1 /. w -> c e]} /. Quiet@FindRoot[(eq1 /. w -> c e) == 0, {e, .1}, Method -> {"Newton", "StepControl" -> "TrustRegion"}], {a, .001, .999, .001}], #[[2]] > 0 && #[[2]] < 1 && #[[3]] < 10^-10 &], {c, .02, .22, .05}]; p2 = ListPlot[Table[sol04[[j]][[All, 1 ;; 2]], {j, 5}], AxesLabel -> {"a", "e"}]; Show[p1, p2]  We also can generate 3D plot using code sol = Select[ Flatten[Quiet@ Table[{a, c, e, Abs[eq1 /. w -> c e]} /. Quiet@FindRoot[(eq1 /. w -> c e) == 0, {e, .9}, Method -> {"Newton", "StepControl" -> "TrustRegion"}], {a, .01, .99, .01}, {c, .01, .5, .01}], 1], #[[3]] > 0 && #[[3]] < 1 && #[[4]] < 10^-10 &]; sol0 = Select[ Flatten[Quiet@ Table[{a, c, e, Abs[eq1 /. w -> c e]} /. Quiet@FindRoot[(eq1 /. w -> c e) == 0, {e, .1}, Method -> {"Newton", "StepControl" -> "TrustRegion"}], {a, .01, .99, .01}, {c, .01, .5, .01}], 1], #[[3]] > 0 && #[[3]] < 1 && #[[4]] < 10^-10 &]; ListPointPlot3D[{sol0[[All, 1 ;; 3]], sol[[All, 1 ;; 3]]}, AxesLabel -> {"a", "c", "e"}, PlotTheme -> "Marketing", PlotRange -> All]  Update 2 The original ppp's code could be modified as follows Clear["Global*"]; L = (w/(b (a e)^b))^(1/(b - 1)); v = -(((-1 + e - i) (-i + (-1 + e + e i) L n) (1/((-1 + e) (-1 + e - i) w) + w/(-1 + e - i) + ((-1 + e) (1 - e L n) (-(1 + i)^2 (-1 + (1 + i)^((1 - e L n)/(L n - e L n)))^2 + a^2 i^2 (-1 + (1 + i)^(1 + 1/(L n - e L n)))^2 W^2))/(a i (1 + i) (1 - e + i) (1 - (1 + i)^((1 - e L n)/(L n - e L n))) (1 - (1 + i)^(1 + 1/(L n - e L n))) (i - (-1 + e + e i) L n) W)))/(i (1 + i + L n + e (-1 + (-2 + e - i) L n)))); z = -((a i (1 + i)^((1 + e)/(1 - e)) (-1 + (1 + i)^((1 - e L n)/(L n - e L n))) (-1 + (1 + i)^(1 + 1/(L n - e L n))) L n W + a (-1 + e) i (1 + i)^((1 + e)/(1 - e)) (-1 + (1 + i)^((1 - e L n)/(L n - e L n))) (-1 + (1 + i)^(1 + 1/(L n - e L n))) L n w^2 W + (-1 + e - i) (1 + i)^(-((2 e)/(-1 + e))) (-1 + e L n) w (1 + 2 i + i^2 - a^2 i^2 W^2 - 2 (1 + i)^(1 + 1/(L n - e L n)) ((1 + i)^(2 + 1/(-1 + e)) - a^2 i^2 W^2) + (1 + i)^(2 + 2/(L n - e L n)) ((1 + i)^((2 e)/(-1 + e)) - a^2 i^2 W^2)))/(a i^2 (1 + i) ((1 + i)^(e/(1 - e)) - (1 + i)^(1/(L n - e L n))) ((1 + i)^(e/(1 - e)) - (1 + i)^((1 + L n)/(L n - e L n))) (1 + i + L n + e (-1 + (-2 + e - i) L n)) w W)); icSOL[x_] := Module[{rule = {n -> 1, W -> 60, i -> 1/10, b -> 6/10, a -> x}}, myEQ = 1/((1 - e)^2 w) - v + z /. {e -> e[w]}; icEQ = D[myEQ, w] /. {e'[w] -> e[w]/w} /. {e[w] -> e}; s = NMinimize[{icEQ^2 /. rule, {0 < e < 1, 0 < w < 2}}, {e, w}]; s] (*Compute initial condition*) icSOL[0.5] (*Out[]= {2.18427*10^-24, {e -> 0.359048, w -> 0.112498}}*) ic = {e, w} /. %[[2]]; eqs = Module[{rule = {n -> 1, W -> 60, i -> 1/10, b -> 6/10}}, eq = 1/((1 - e)^2 w) - v + z; rhs = D[{eq, {e, w} . D[eq, {{e, w}}]} /. {w -> w[a], e -> e[a]}, a]; Thread[rhs == {0, 0}] /. rule]; sol = NDSolveValue[{eqs, e[1/2] == ic[[1]], w[1/2] == ic[[2]]}, {e, w}, {a,$MachineEpsilon, 1},
Method -> {"EquationSimplification" -> "Residual"}]


Note that we have solution for 0<a<0.504197417091105 only. Visualization

{Plot[sol[[1]][a], {a, 0, .5}, PlotRange -> {0, 1},
AxesLabel -> {a, "e"}],
Plot[sol[[2]][a], {a, 0, .5}, PlotRange -> All,
AxesLabel -> {a, "w"}],
Plot[sol[[1]][a]/sol[[2]][a], {a, 0, .5}, AxesLabel -> {"a", "e/w"}],
Plot[a  sol[[1]][a]/sol[[2]][a], {a, 0, .5},
AxesLabel -> {"a", "a e/w"}]}


• Thanks! Some questions: 1) It seems s2 is a typo? 2) It seems c is considered as a constant parameter just like b, i, W, n. But it is not. It is endogenously determined just like e and w. Maybe I'm misreading your code?
– ppp
Commented May 31 at 18:15
• So for point (2) in my above comments, here is a clarification of how c is endogenous in my system: From the first equation in my post, e=e(w; a,i,b,n,W) can be derived, which turns out to be a concave and increasing function in w. Now the (e,w) pair that satisfies the second equation in my post would occur where a linear ray from the origin touches the e=e(w) curve. And c would refer to the slope of the ray. In this context, an increase in a will shift the e=e(w) curve to the right and hence will make the tangent line flatter, i.e. making c smaller.
– ppp
Commented May 31 at 18:23
• You can take a look at the first two graphs in the link below and you'll see what I mean. mathematica.stackexchange.com/questions/287287/…
– ppp
Commented May 31 at 18:28
• @ppp 1) First code generates optimal solution only. 2) Second code generates solution for different c, and last picture shows this kind of solution. 3) If you suppose that $c=c(a)$ then you need additional equation to compute this function. Another words you need initial condition for differential equation. Commented Jun 1 at 4:52
• @ppp Please, see Figure 1 in my answer above. Solution has two branches connecting in critical point. In this regards it could be better to use e as a variable and a[e],w[e] as functions. Commented Jun 5 at 2:37