2
$\begingroup$

Consider the following non-linear equation:

$$ k x = \frac{\frac{9}{4}-y^{22/25}}{2\ 2^{19/25} y^{22/25}-\frac{9}{4} \left(10 y+e^{-y}\right)^{22/25}-1} $$

where k is a constant, i.e, $k=0.3516$, $y=y(x)$ and $y(0)=\log(2.25^{0.88}) \quad x \in [0,1]$

I want to solve the numerical solution of $y$.

My trail as shown below:

$$ \frac{d(k x)}{dx} = \frac{d}{dx}\left(\frac{\frac{9}{4}-y^{22/25}}{2\ 2^{19/25} y^{22/25}-\frac{9}{4} \left(10 y+e^{-y}\right)^{22/25}-1}\right) $$

rightHand = 
  D[(-y[x]^(22/25) + 9/4)/(-9/4*(10 y[x]/1 + E^(-y[x]))^(22/25) + (4*y[x])^(22/25) - 1), x] // Simplify;

leftHand = 0.3516;
sol = NDSolve[{leftHand == rightHand, y[0] == Log[2.25^0.88]}, y, {x, 0, 1}]

However, Mathematica gives the following error:

NDSolve::ndsz: At x == 0.5671511305713252`, step size is effectively zero; singularity or stiff system suspected. >>

Plot[Evaluate[y[x] /. sol], {x, 0, 1}, PlotRange -> All]

enter image description here

$\endgroup$
13
  • $\begingroup$ I found a solution by writting the equation as y'[x]==newRightHandSide using eq = First@First@Solve[leftHand == rightHand, Derivative[1][y][x]] /. Rule -> Equal; $\endgroup$ Dec 8, 2022 at 17:25
  • $\begingroup$ then sol = NDSolve[{eq, y[0] == Log[2.25^0.88]}, y, {x, 0, 1}] but the solution has a strong slope as shown in this plot . Is that expected ? $\endgroup$ Dec 8, 2022 at 17:33
  • $\begingroup$ Also having used this method in systems of non linear equations this stiff message can occur at times. For the case of multiple non linear equations I usually used one of the options of NDSolve for stiff problems but I do not remember the name now. You can maybe check the documentation for stiff problems if that happens in the future. $\endgroup$ Dec 8, 2022 at 18:07
  • $\begingroup$ @userrandrand Thanks, I notice that your solution still has the error NDSolve::ndsz: At x == 0.5671511305713252, step size is effectively zero; singularity or stiff system suspected. >> $\endgroup$
    – Lavender
    Dec 8, 2022 at 22:57
  • $\begingroup$ I did not receive any error $\endgroup$ Dec 8, 2022 at 22:58

1 Answer 1

3
$\begingroup$

Try this:

f[y_] := (-y^(22/25) + 
    9/4)/(-9/4*(10 y/1 + E^(-y))^(22/25) + (4*y)^(22/25) - 1);

lst = Table[{kx, 
    FindRoot[f[y] == kx, {y, 0.1 + 2 kx}][[1, 2]]}, {kx, -0.07, 0.07, 
    0.001}] // Chop;

yielding

ListLinePlot[lst, 
 AxesLabel -> {Style["y", 16, Italic], Style["kx", 16, Italic]}]

enter image description here

Edit: There is a more direct way of resolving the plot by using a ParametricPlot:

ParametricPlot[{(-y^(22/25) + 
     9/4)/(-9/4*(10 y/1 + E^(-y))^(22/25) + (4*y)^(22/25) - 1), 
  y}, {y, 0, 5}, AspectRatio -> 0.7, 
 AxesLabel -> {Style["y", 16, Italic], Style["kx", 16, Italic]}]

enter image description here

In this sense, one can regard your equation as already resolved with respect to kx in the parametric form.

Have fun!

$\endgroup$
3
  • $\begingroup$ Thanks for your solution. I have a question with {kx, -0.07, 0.07, 0.001}. Because $x \in [0,1]$, $k x \in [0, k]$. $\endgroup$
    – Lavender
    Dec 8, 2022 at 22:55
  • $\begingroup$ from the data in the ListPlot that is y as a function of kx no ? @Lavender looking at the plot of kx as a function of y , one sees a plateau for y large which seems to infer that y as a function of kx, the reciprocal, diverges.Specifically, it seems to diverge around 0.07 which might be why that was chosen above. $\endgroup$ Dec 9, 2022 at 0:12
  • $\begingroup$ @Lavender You just take those values that you are looking for. That's up to you. There is a more direct way of doing this. See the later edit. $\endgroup$ Dec 9, 2022 at 10:27

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.