calculation variable x besselj

3 views (last 30 days)
Thomas Krone
Thomas Krone on 17 Oct 2017
Edited: John Kelly on 18 Oct 2017
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?

Accepted Answer

KSSV
KSSV on 17 Oct 2017
eqn = @(x) (x.*besselj(1, x))./besselj(0, x)-5 ;
x0 = [1 4 7 10 12 15] ;
x = fsolve(eqn,x0)
  3 Comments
John D'Errico
John D'Errico on 17 Oct 2017
Edited: John D'Errico on 17 Oct 2017
@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.

Sign in to comment.

More Answers (3)

John D'Errico
John D'Errico on 17 Oct 2017
Edited: John D'Errico on 18 Oct 2017
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.
  1 Comment
KSSV
KSSV on 18 Oct 2017
+1 for your in depth solution....

Sign in to comment.


Torsten
Torsten on 17 Oct 2017
and put the mouse over the red points in the plot.
Enough solutions ?
Best wishes
Torsten.

John BG
John BG on 17 Oct 2017
Edited: John Kelly on 18 Oct 2017
Hi Thomas
One way to solve this equation is with the function intersections.m by Mr Schwarz's available here
.
.
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
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);
  3 Comments
John D'Errico
John D'Errico on 18 Oct 2017
Edited: John Kelly on 18 Oct 2017
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.

Sign in to comment.

Products

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!