calculation variable x besselj

7 Ansichten (letzte 30 Tage)
Thomas Krone
Thomas Krone am 17 Okt. 2017
Bearbeitet: John Kelly am 18 Okt. 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?

Akzeptierte Antwort

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

Melden Sie sich an, um zu kommentieren.

Weitere Antworten (3)

John D'Errico
John D'Errico am 17 Okt. 2017
Bearbeitet: John D'Errico am 18 Okt. 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;
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.
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?
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 =
I've started fzero out in a bad place. Why do I know that?
ans =
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 =
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 =
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 =
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;
Did it work? Yes. It worked because I was careful.
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.
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 =
ans =
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 =
ans =
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 Kommentar
KSSV am 18 Okt. 2017
+1 for your in depth solution....

Melden Sie sich an, um zu kommentieren.

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

John BG
John BG am 17 Okt. 2017
Bearbeitet: John Kelly am 18 Okt. 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;
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 =
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
  3 Kommentare
John D'Errico
John D'Errico am 18 Okt. 2017
Bearbeitet: John Kelly am 18 Okt. 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.

Melden Sie sich an, um zu kommentieren.


Community Treasure Hunt

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

Start Hunting!

Translated by