Non-linear equation of functions

Question

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);
areasm = Pi radiussm^2; 
dareasm = D[Pi radiussm^2, z];

DIFFERENTIAL EQUATION:

rsol = NDSolve[{((k1^2)/((areasm^2) r[z]) - (k2 r[z]^gamma) gamma) r'[z] == -((k1^2)/(areasm^3)) dareasm, r[0] == 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!

Answer 1

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][[2]]], {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.

Answer 2

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.

Answer 3

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}]