Solving non-linear equation with besselj

Hello,

I am trying to solve this equation below and hope to get a closed form so I can plot it. I understand that due to the periodic nature of Bessel functions, "solve" does not work and I can actually get results using vpasolve with good guesses; or simply use a graphic method manually to make sure the solution coincides with the first root of the Bessel. However, I suspect there is a way to get symbolic solution using solve so it spits out a closed form rather than a number if I only restrict solutions to the first zero of the besselj I use. Is this possible?

Thank you for your input!

syms f 
assume (f > 1e6)
r=20e-6;
n=2;
E=28e9; 
rho=8095; 
sigma=0.22; 

w=2*pi*f;
zeta=r*sqrt(2*rho*w.^2.*(1+sigma)/E);
q=zeta.^2./(2*n^2-2);
xi=sqrt(2/(1-sigma));

eqn=(Psi_2(zeta/xi)-2-q )*(2*Psi_2(zeta)-2-q) == (n*q-n)^2 ; %Equation I want to solve with only 1st root of Bessel function
f=1e-6*solve(eqn,f);

% Definition of Psi_2 function 
function f = Psi_2(theta) 
    f=theta.*besselj(1,theta)./besselj(2,theta);
end

Answer 1

[Walter Roberson]

I would very much doubt you can find a symbolic solution. You are trying to solve

49379500000000 * besselj(2, 1/35000000000*98759^(1/2)*14^(1/2)*Z*pi) * ...
besselJ(1,1/350000000000*39^(1/2)*98759^(1/2)*14^(1/2)*Z*pi) * ...
39^(1/2)*98759^(1/2)*14^(1/2)*pi^2*Z^2 + ...
987590000000000 * besselj(1, 1/35000000000*98759^(1/2)*14^(1/2)*Z*pi) * ...
besselj(2, 1/350000000000*39^(1/2)*98759^(1/2)*14^(1/2)*Z*pi) * ...
98759^(1/2)*14^(1/2)*pi^2*Z^2 + ...
9753340081 * besselj(2, 1/35000000000*98759^(1/2)*14^(1/2)*Z*pi) * ...
besselj(2,1/350000000000*39^(1/2)*98759^(1/2)*14^(1/2)*Z*pi) * pi^3*Z^3 - ...
20739390000000000000000000 * besselj(1, 1/35000000000*98759^(1/2)*14^(1/2)*Z*pi) * ...
besselj(1, 1/350000000000*39^(1/2)*98759^(1/2)*14^(1/2)*Z*pi) * 39^(1/2)*pi*Z + ...
525000000000000000000000000000 * besselj(2, 1/35000000000*98759^(1/2)*14^(1/2)*Z*pi) * ...
besselj(1, 1/350000000000*39^(1/2)*98759^(1/2)*14^(1/2)*Z*pi) * 39^(1/2)*98759^(1/2)*14^(1/2) + ...
10500000000000000000000000000000 * besselj(1, 1/35000000000*98759^(1/2)*14^(1/2)*Z*pi) * ...
besselj(2, 1/350000000000*39^(1/2)*98759^(1/2)*14^(1/2)*Z*pi) * 98759^(1/2)*14^(1/2) - ...
207393900000000000000000000 * besselj(2, 1/35000000000*98759^(1/2)*14^(1/2)*Z*pi) * ...
besselj(2, 1/350000000000*39^(1/2)*98759^(1/2)*14^(1/2)*Z*pi)*pi*Z

for the set of values of Z such that the expression becomes 0.

That is not just a single zero of a bessel function: there are 12 BesselJ functions there. And there is no symbolic solution to zeros of even a single BesselJ .

[Ozan Erturk]

I wanted to get an expression for f vs. r in the script I just let r=20e-6. I thought it would be possible if there is a way to force each besselj function to restict to it's first zero only, but I am guessing now that it is possible, right?

I tried vpasolve with an intial guess and it gives me the expected answer. However, sometimes depending on the initial guess, I don't get the first solution even though it is still numerically close to the first solution. Is it normal for vpasolve to output in this way? If so, then I want to approach it graphically to make sure I am extracting the correct solution, not the higher order ones.

Here is my second approach:

If I plot/solve the LHS and RHS of the equation separately, the solution I want is the first crossing only. I do not care about the higher order solutions. So if I numerically find the first point where LHS=RHS I am done (when I repeat this for the range of "r" values I want). Does this make sense?

Thank you for your time Walter!

-Ozan

[Walter Roberson]

Once you have solved one, pass that result in as the initial guess to the next one.

As I test with those equations, I find that the number of solutions is sensitive to the value of r and f.

[Ozan Erturk]

Thanks! That seems to be the only way I guess, at least the variation of f can be guessed from the physics of it as well.

Thank you for your input!

[Walter Roberson]

The below has to do with the solution. If you take lhs(eqn)-rhs(eqn) and solve() that for f, using a better symbolic package, then you get a RootOf() for quite similar to what I posted first -- a function of _Z such that the function must equal 0 at _Z, and _Z is then the frequency, f, at which the equation holds true.

If you remove the RootOf() level and plot just the internal part, the part that needs to become 0, then the places where it does become 0 are the solutions. 0 is the flat horizontal plane in the image, so there are solutions everwhere the curve intersects the flat plane.

So the smaller the r, the higher the lowest frequency of solution. With high enough r, the solutions start coming in pretty close (e.g., as r approaches 1, solutions down to f = 600-ish appear.)

This is to be expected: the besselJ calls end up having f*r in them, so if you were to create a new variable FR = f*r then you get to a 1D equation. If you know a solution for one f*r value then you can predict the solution for a different r by inverse proportions.

[Ozan Erturk]

Although I am not sure, I think I understand the problem, and your suggestion. I implemented it and called a variable "zeta" which has both f and r. This was the original function without any f and r substitutions in fact. However, I don't think solve can still handle this. Do I need to modify "solve" ? or any other symbolic solver in matlab that I don't know about maybe?

syms zeta positive

n=2;
sigma=0.22; %Poisson's ratio
q=zeta^2/(2*n^2-2);
xi=sqrt(2/(1-sigma));
eqn=(Psi_2(zeta/xi)-2-q )*(2*Psi_2(zeta)-2-q) == (n*q-n)^2 ;
zeta=solve(eqn,zeta)

Warning: Unable to find explicit solution. For options, see help.

 
zeta =
 
Empty sym: 0-by-1
 

function f = Psi_2(theta)
    f=theta.*besselj(1,theta)./besselj(2,theta);
end