# Non-linear equation of functions

I would like to solve (I think with NSolve, but I am not sure) and plot the function $$w(z)$$: the relation contains other known functions dependent on $$z$$, which are areasm and p, where p is the function resulting from an NDSolve. In this relation there is no an explicit dependence of $$w$$ from $$z$$.

I try to write the wall code for clarity:

PARAMETERS:

d1 = 0.8*^-3;
d2 = 0.5*^-3;
l = 17*^-3;
ltot = 22*^-3;
l1 = 13.8*^-3;
l2 = 14.5*^-3;
aperture = 0.3;
slope = 10;
p1 = 2;
T = 300;
vi=1*^3;
step = 10^3;
Kb = 1.380649 10^(-23);
gamma = 1.660;
R = 8.314462618;
vsound = 965;
area1 = Pi (d1/2)^2;
pi = p1*10^6;
ri = pi/(Kb T);
k1 = ri vi area1;
k2 = pi ri^(-gamma);


FUNCTIONS:

radiussm = ((((d1/2) Exp[l1 slope*10^3])+((d2/2) Exp[slope*10^3 z]))/(Exp[l1 slope*10^3]+Exp[slope*10^3 z]))+UnitStep[z-l2] aperture (z-l2);


DIFFERENTIAL EQUATION:

rsol = NDSolve[{((k1^2)/((areasm^2) r[z]) - (k2 r[z]^gamma) gamma) r'[z] == -((k1^2)/(areasm^3)) dareasm, r == ri},r,{z, 0, ltot},InterpolationOrder -> All];
p = k2 (r[z]^gamma);


PROBLEM:

velocity = NSolve[k1 == ((areasm*(k2*((First[Evaluate[p[z] /. rsol]])^gamma)))/Sqrt[T])*(Sqrt[gamma/R])*(w/vsound)*(1 + ((gamma - 1)/2)*(w/vsound)^2)^(-(gamma + 1)/(2*gamma - 2)), w];
velocitytab = Table[w /. velocity, {z, 0, ltot, ltot/step}]
Plot[Evaluate[w /. velocitytab], {z, 0, ltot}, PlotRange -> All,AxesLabel -> {"z [m]", "v [m/s]"}, ImageSize -> 350]


...but I didn't get any plot! Where did I go wrong? Thank you in advance!

• How long does velocity = NSolve[...] run? The calculation runs quite long on my machine. – Max1 May 1 at 18:15
• Yes, it takes a really long time because of the two gamma-dependent exponents. If I enter integers instead of them, the NSolve calculation takes a few seconds, but it still has other problems, at least in the plot. However I do not understand why the gamma in the NSolve creates all this slowdown, because it appears as an exponent also in the previous NDSolve which has no problems instead. – gecos May 1 at 20:02
• I don't quite understand the argument of NSolve. It depends on z, but its not a numeric value when z is passed as an argument. A second thought: You should be able to do the NSolve calculation inside Table for each z value right? Try to calculate that for like 10 value and then look if you actually get numeric values or garbage. – Max1 May 1 at 20:58
• In the NSolve argument: I would like to find the w(z) function. It depends on z as it depends on the areasm and p function that both are functions of z. z for me are simply steps of an interval between z=0 and z=ltot, as it is in the NDSolve argument. However, I was not able to enter this range of z values in the NSolve. – gecos May 1 at 21:35
• I cannot calculate the NSolve either outside or inside the Table I believe because of the fractional exponents. If I insert some integer coefficients, the Table gives the same result both if I insert the NSolve in the argument of the Table, and if I only recall the result of the NSolve in the Table. I can't understand if they are numerical values or not from the Table result, well I think they are, but then I don't understand why the plot is empty. – gecos May 1 at 21:44

Take the definitions from the question up until the section PROBLEM:

First have a look at the definition of velocity=NSolve[...]. The argument of NSolve is a function if w and depends on z. A range of z values is given (0 to ltot). Let's pick one of the values and perfrom the root search for a fixed z. I suppose the argument of the NSolve can then be written as

velocityRoots[z0_] = (((areasm*(k2*((First[p /. rsol])^gamma)))/Sqrt[T])*(Sqrt[gamma/R])*(w/vsound)*(1 + ((gamma - 1)/2)*(w/vsound)^2)^(-(gamma + 1)/(2*gamma - 2)) - k1) /. z -> z0;


Note that I have replaced First[Evaluate[p[z]/.rsol]] (from the original question) with First[p /. rsol], as I believe that is what is actually intended. Now the upper function is the function which we wish to know the roots of. Let's inspect that function. velocityRoots[ltot/100] gives

$$-2.42714\times 10^{23}+\frac{3.92562\times10^{-39}w}{\left(1+3.54372\times10^{-7} w^2\right)^{2.01515}}$$

By inspection the magnitude of the offset term is huge and the magnitude of the function which follows is tiny.

To convince graphically let's plot the function:

Plot[Evaluate[velocityRoots[ltot/100][]], {w, -10000, 10000},PlotPoints -> 200] Due to offset of magnitude $$10^{23}$$ there will be no roots for the given z value. Unfortunately if you play around with other z values between 0 and ltot the situation does not change. The offset is huge and the magnitude of the function is tiny. The overall function is almost constant.

In conclusion: Mathematica cannot plot anything because your function does not have any roots at all.

• Ok, I understood what the problem was, in fact the code didn't work because it was really badly defined in physical problem because of errors in the definition of the function to solve in the NSolve. The fact of plotting the roots for a z0 helped me a lot! thank you! – gecos May 5 at 11:18
• @gecos your welcome, as a general advice: It is often difficult to get all parameters right when modeling a complex physical system. Dividing the problem in smaller chunks, which can be checked more easily, can be immensely helpful. Personally I have taken a few habits from coding like writing a small design doc and creating small mile stones even for simulation projects. Visualizing work can also be very helpful for oneself. – Max1 May 5 at 11:28

PROBLEM can be solved the way put in the question.

Since this is numerical problem get rid of all these symbolic constants and use numbers instead.

p[z] and Evaluate are redundant.

There might be a solution existing on {w,z} but this is hard for Mathematicas methods to find.

Informativ is the message if Solve is used instead of NSolve: The problem is numerical exactness. Mathematica is famous for exact number representation but that is a hard problem for numerics.

So make up a mind about precision and continue with: The Solve::inex message has no page for itself.

It might be of use to restrict the evaluation in z on the domain of the solution rsol.

Since z is a general plugin the rsol into the solution with numerical Precision. A trick in general is to used Rationals as constants. These have the desired infinite precision. Mathematica will internally optimize the results for that input.

Again, in general, avoid options for NDSolve like InterpolationOrder or PrecisionGoal until at least there is a proper error estimate or calculation or the differential equation is safe to the solved with these options. There ongoing discussions on whether the effects are discussed in the Mathematica documentation or not. Since everything is somehow hidden or inherent and error are often hidden there is a long way path to understand my advice from the Mathematica documentation.

Not everybody knows but those who search more in the documentation know that there is NDSolve and NDSolveValue and ParametricNDSolve and some more built-ins for solve specials paths for solutions.

Thanks to the answers I received, now the code works and does what it has to do, that is, the NSolve finds roots even for fractional exponents: I leave here a simplified example of the use of NSolve to find the roots of functions dependent on other functions of the same independent variable, and the related plots. Suppose we have a parabolic function with a coefficient a, which depends on a variable y. We plot the roots of the parabola for a given y, and then the roots with respect to y

a = y^2
d1 = x^2 + a 2 x - 3;
Dp = First[NSolve[d1 == 0, x]]
dd = d1 /. y -> 2
Plot[dd, {x, -10, 10}]
Plot[Evaluate[x /. Dp], {y, 1, 4}]