calculation variable x besselj

For heat transfer calculations where I'm working on I have to solve the equation:

(x*besselj(1, x))/besselj(0, x) = 5

I first tried it using WolframAlpha,where I got 4 values of x. (Namely, x=1.989; x=4.713; x=7.618; x=10.622)

Because I would like to have 6 values of x, i will use Matlab.

Is there anybody that can tell me how to solve the equation using Matlab?

Answer 1

[KSSV]

eqn = @(x) (x.*besselj(1, x))./besselj(0, x)-5 ;
x0 = [1 4 7 10 12 15] ;
x = fsolve(eqn,x0)

[Thomas Krone]

Thanks for your answer!

Is it also possible to make '5' variable like

z = (0:0.2:1)';
eqn = @(x) (x.*besselj(1, x))./besselj(0, x)- (z)  ;
x0 = [1 4 7 10 12 15] ;
x = fsolve(eqn,x0)

When I do this there is no solution found.

[Stephen23]

@Thomas Krone: you could put it into a loop:

z = 0:0.2:1;
x0 = [1,4,7,10,12,15];
m = z(:)*x0 + NaN; % preallocate output
for k = 1:numel(z)
  eqn = @(x) (x.*besselj(1, x))./besselj(0, x)- z(k);
  m(k,:) = fsolve(eqn,x0);
end

But you should definitely read John D'Errico's answer too!

[John D'Errico]

@Thomas - If you want to solve for a variety of targets you cannot set z to a vector there, at least not without writing code to do it properly, as I showed.

You need to recognize that for ANY value of z as a target, there will be infinitely many solutions. So if you wanted to solve for the first 10 solutions for each value of z=(0:0.2:1)', the result must in fact be a 6x10 array of solutions.

I do show how to solve the problem for various targets in my answer.

Answer 2

[John D'Errico]

Note that there are infinitely many solutions to this problem. You can see that if you examine the plot shown.

ytarget = 5;
bfun = @(x,ytarget) (x.*besselj(1, x))./besselj(0, x) - ytarget;
ezplot(@(x) bfun(x,ytarget),[0,100])
grid on

In that plot, you are looking for a zero crossing, since I subtracted the target value.

I wrote bfun this way so that you can see you can change the target value for the value of the function easily. You did ask for that in one of your comments. Also note that I used the .* and ./ operators to make evaluation of this function vectorized.

Can you find as many solutions as you wish? Yes. That is pretty easy in fact, since you should recognize that the solutions will be separated by a near constant increment, and that increment approaches pi from below. You need to be careful, since fzolve or fzero can become confused. Why?

This function has singularities wherever besselj(0,x) is zero. In fact, every true solution will be separated by such a singularity.

So, lets find the first 10 true solutions.

xroots = zeros(1,10);
x0 = 2;
for i = 1:10
  xroots(i) = fzero(@(x) bfun(x,ytarget),x0);
  x0 = xroots(i) + 2.5;
end

xroots
xroots =
  1.9898  4.7131  7.6177  10.622  13.679  16.763  19.864  22.975  26.094  29.217

Are they true solutions? Yes.

bfun(xroots,ytarget)
bfun(xroots,ytarget)
  ans =
  -3.5527e-15   4.4409e-15  -1.0658e-14   2.4869e-14   3.6415e-14  -5.1514e-14  -2.8422e-14   3.0198e-14  -3.3751e-14    3.908e-14

What is the difference between consecutive roots as found?

diff(xroots)
ans =
  2.7233  2.9046  3.0046  3.0563  3.0844  3.101  3.1114  3.1183  3.1231

This is why I incremented x0 by only 2.5 in this line, even though I know that the increment is closer to pi. By incrementing x0 by a number that sufficiently is smaller than pi, but also not too small, I have forced fzero to always look in a spot that will allow it to converge to the next root of interest.

    x0 = xroots(i) + 2.5;

Had I chosen a too large increment there, fzero would have found some spurious roots, where the function crosses a singularity. For example, suppose I do this?

format long g
ytarget = 5;
xbad = fzero(@(x) bfun(x,ytarget),2.5)
xbad =
        2.40482555769577

I've started fzero out in a bad place. Why do I know that?

bfun(xbad,ytarget)
ans =
      -3.4020364853122e+15

So this is not a solution to the problem. fzero found a zero crossing, but it found a crossing where the function passes through a singularity. The denominator of your function crosses zero near that location. So the function transitions from inf to -inf, crossing zero but never passing through the value we need. Fzero got confused.

syms x, vpasolve(besselj(0,x),2)
ans =
2.4048255576957727686216318793265

As it turns out, fsolve will not fall into that trap, because it is a different class of solver than fzero.

xgood = fsolve(@(x) bfun(x,5),2.5,opts)
xgood =
          4.71314228694607

Addendum:

Thomas wants to find a solution for various target values. Can this be done in a vectorized form? Thus...

ytarget = (0:0.2:1)';

The problem is there are 6 values here. For each target value, we can find infinitely many roots.

ytarget = (0:.2:1)';
bfun = @(x,y0) (x.*besselj(1, x))./besselj(0, x) - ytarget;
fsolve(@(x) bfun(x,ytarget),repmat(.2,6,1))

Equation solved.

fsolve completed because the vector of function values is near zero
as measured by the default value of the function tolerance, and
the problem appears regular as measured by the gradient.

<stopping criteria details>

ans =
    0.0058484
      0.61697
      0.85158
       1.0184
        1.149
       1.2558

But that is only the first solution. Now, you can carefully choose an increment for fsolve to find the next root.

ytarget = (0:.2:1)';
bfun = @(x,y0) (x.*besselj(1, x))./besselj(0, x) - ytarget;
nroots = 10;
xroots = zeros(numel(ytarget),nroots);

opts = optimset('fsolve');
opts.Display = 'none';

x0 = repmat(.5,size(ytarget));
for i = 1:nroots
  xroots(:,i) = fsolve(@(x) bfun(x,ytarget),x0,opts);
  % increment by a bit less than pi each time.
  x0 = xroots(:,i) + pi - .25;
end

Did it work? Yes. It worked because I was careful.

xroots
xroots =
     0.006345       3.8317       7.0156       10.173       13.324       16.471       19.616        22.76       25.904       29.047
      0.61697       3.8835        7.044       10.193       13.339       16.483       19.626       22.769       25.911       29.054
      0.85158       3.9344       7.0723       10.213       13.354       16.495       19.636       22.778       25.919       29.061
       1.0184       3.9841       7.1004       10.232       13.369       16.507       19.646       22.786       25.927       29.067
        1.149       4.0325       7.1282       10.252       13.384       16.519       19.657       22.795       25.935       29.074
       1.2558       4.0795       7.1558       10.271       13.398       16.531       19.667       22.804       25.942       29.081
>> bfun(xroots,ytarget)
ans =
   2.0129e-05   5.2795e-08   3.2883e-11   7.3099e-12   4.0079e-12   2.2535e-12   1.7552e-12   1.3051e-12   8.1562e-13   6.2715e-13
   2.7756e-17   3.3307e-16   9.1267e-12   7.3611e-12    5.862e-12   4.7589e-12   4.0418e-12   2.6944e-12   2.4602e-12   1.9923e-12
            0  -4.4409e-16   5.0898e-13   4.0423e-12   6.0018e-12   5.7701e-12   5.4801e-12   4.7149e-12   4.2417e-12   3.8285e-12
            0   1.5543e-15   7.7716e-16   1.3242e-12   4.4678e-12   6.1195e-12   6.4949e-12   6.6159e-12    5.992e-12   5.4791e-12
  -1.1102e-16  -5.5511e-16   7.8737e-13   1.6864e-13   2.6555e-12   5.5311e-12   6.9458e-12   7.7957e-12   7.7064e-12   7.6598e-12
            0  -8.8818e-16   5.8382e-12  -5.9952e-15   9.7367e-13   3.7421e-12   6.6613e-12    8.257e-12   9.1596e-12   9.0155e-12

As you can see, the difference between consecutive roots is again roughly pi as the root index increases.

diff(xroots,[],2)
ans =
     3.8254       3.1839       3.1579       3.1502       3.1469       3.1452       3.1442       3.1436       3.1432
     3.2665       3.1605       3.1491       3.1456       3.1441       3.1433       3.1428       3.1425       3.1423
     3.0828        3.138       3.1404        3.141       3.1412       3.1413       3.1414       3.1415       3.1415
     2.9656       3.1163       3.1318       3.1364       3.1384       3.1394         3.14       3.1404       3.1407
     2.8835       3.0957       3.1234       3.1319       3.1356       3.1375       3.1386       3.1393       3.1398
     2.8237       3.0763       3.1152       3.1274       3.1328       3.1356       3.1372       3.1383        3.139

Finally, expect some inaccuracy in a root finder when the target on this problem is at zero.

ezplot(@(x) bfun(x,0),[-.1 .1])

So here, the function has what would be properly called a higher order root. Here the root is of multiplicity 2, since both the function and the derivative are zero there. This causes the solution to be poor. Worse, tools like fzero may never actually find that root, even though 0 is a solution.

fzero(@(x) bfun(x,0),.1)
ans =
    -2.4048

bfun(0,0)
ans =
     0

There is no zero crossing at that point, and fzero looks for a zero crossing. fsolve can succeed, since it is essentially a Newton-like scheme.

fsolve(@(x) bfun(x,0),.1,opts)
ans =
     0.012521

bfun(ans,0)
ans =
   7.8384e-05

But it was not very close to the true solution at zero, again, because of the quadratic nature of the problem in that vicinity.

[KSSV]

+1 for your in depth solution....

Answer 3

[Torsten]

http://www.wolframalpha.com/input/?i=Solve%5Bx*J1(x)%2FJ0(x)+%3D%3D+5,x,Reals%5D

and put the mouse over the red points in the plot.

Enough solutions ?

Best wishes

Torsten.

Answer 4

[John BG]

Hi Thomas

this is John BG <mailto:[e-mail] [e-mail]>

One way to solve this equation is with the function intersections.m by Mr Schwarz's available here

.

https://uk.mathworks.com/matlabcentral/fileexchange/11837-fast-and-robust-curve-intersections?s_tid=srchtitle

.

intersections.m accurately gets all zeros within the range you specify, as long as one doesn't input too many points:

clear all;clc;close all;
L=30;dx=.005;x=[-L:dx:L];
y=(x.*besselj(1, x))./besselj(0, x) - 5;
plot(x,y);grid on

.

.

zoom detail of the zeros

.

.

xc=intersections(x,y,[-L L],[0 0])

xc =
 -29.216809655219919
 -27.491519805992304
 -26.093677068490667
 -24.352527057988961
 -22.975360511598623
 -21.213361945994276
 -19.863965163930313
 -18.073934761523411
 -16.762981766775486
 -14.934080915449329
 -13.678556147466056
 -11.793463084047159
 -10.622297193226389
  -8.651269069199444
  -7.617703469680007
  -5.524921507288155
  -4.713136132229037
  -2.400172515733278
  -1.989812636519481
   1.989812636519481
   2.400172515733278
   4.713136132229037
   5.524921507288155
   7.617703469680007
   8.651269069199444
  10.622297193226389
  11.793463084047159
  13.678556147466056
  14.934080915449329
  16.762981766775486
  18.073934761523411
  19.863965163930313
  21.213361945994276
  22.975360511598623
  24.352527057988961
  26.093677068490667
  27.491519805992304
  29.216809655219919

.

.

Mr Krone, if you find this answer useful would you please be so kind to consider marking my answer as Accepted Answer?

To any other reader, if you find this answer useful please consider clicking on the thumbs-up vote link

thanks in advance

John BG

<mailto:[e-mail] [e-mail]>

Additional comment:

Also tried the following but the following lines of work do not get the accuracy of function intersections.m

find(y(y==0))

find(y(abs(y)<0.00001))

N1=find(y(y>-.001));N2=find(y(y<.001));intersect(N1,N2)
N1=find(y(y>-.0001));N2=find(y(y<.0001));intersect(N1,N2);
N1=find(y(y>-.00001));N2=find(y(y<.00001));intersect(N1,N2);
N1=find(y(y>-.0000001));N2=find(y(y<.0000001));intersect(N1,N2);
N1=find(y(y>-.000000001));N2=find(y(y<.000000001));intersect(N1,N2);

[John D'Errico]

Yes, it looks vaguely like there is a solution at that location. But is it a solution?

x2 = 2.400172515733278;
bfun = @(x,y0) (x.*besselj(1, x))./besselj(0, x) - 5;
bfun(x2 + [-0.1 0 0.1])
ans =
       17.395       511.32      -30.639

Is there a solution at x2? NO. This is the approximate location of a singularity. The function will never attain the value of 5 in that vicinity.

The problem is that besselj(0,x) has a root at that location.

syms x, vpasolve(besselj(0,x),2)
ans =
2.4048255576957727686216318793265

So the function in question, since it has besselj(0,x) in the denominator will cross zero.

bfun(2.4048255576957727686216318793265)
ans =
   Inf

bfun(2.4048255576957727686216318793265 + 10*eps)
ans =
        -992260641549379

[John D'Errico]

The power of intersections is not in question. It is not perfect. Both intersections AND fzero will return a solution at those singularities. Regardless they are incorrect. A sign change is NOT a root, not when a singularity from -inf to inf occurs at the sign change. The denominator is known to cross at those points. That will cause exactly the class of singularity I have described.

For example, is there a root of tan(x)==0 at x==pi/2? Of course not!

The tangent function shows a sign change there, but there are NO roots. If you would argue that x==pi/2 is a solution to tan(x)==0, then you you could as easily argue that x==pi/2 is a solution to the problem tan(x)==k, for ANY arbitrary value of k. Again, both intersections and fzero will agree with you. But they are clearly wrong in that case.

So, while the use of intersections here is an interesting idea, it is a flawed one.

[John D'Errico]

The use of intersections is fine. In fact, it is a nice alternative way to solve the problem, with the caveat that the solutions produced, even the valid solutions will be moderately poor. The reason is that intersections uses linear interpolation to locate the solutions.

For example, bfun([1.989812636519481,4.713136132229037,7.617703469680007],5) ans = -3.02455609411112e-05 -6.16542254521235e-05 -4.61635674353644e-05

So, while they are approximate, they are only first order approximations, producing a zero with error on the order of 1e-5. These would be great starting values for an iterative root finder.

But anyone who reads this needs to see that intersections produces spurious "solutions", that are not in fact solutions. Just because a discontinuous function crosses zero does not mean that it ever takes on the value zero.