Solve system for Roots of Bessel Function

I have this Bessel function that I am trying to solve for the roots an and bn.They are dependent upon T ,K, B, which I was able to figure out, but I need to be able to use this in my infinite summation model, so being able to solve for more would be useful.

syms a_n b_n T K B

I tried to use vpasolve and fzero function with (T=10,B=1,K=0.5) but didnt work.

What function should i use for solving this system?

Any help would be greatly appreaciated.

Thank you

[Torsten]

What infinite summation model do you refer to in which you need to solve a system of equations containing Bessel functions ?

[mery]

Hi, I need the roos a_n and b_n to plot this sum for n=1000

[Torsten]

In the infinite sum, only the an appear.

Could you give a literature link to this sum representation ?

Answer 1

[Torsten]

This is a very empirical code for your problem. Success is not guaranteed in all cases.

T = 10;
B = 1;
K = 0.5;
fun = @(a,b)[(besselj(0,K*a)+b*bessely(0,K*a))*(1-T*a)+B*a*(besselj(1,K*a)+b*bessely(1,K*a));...
             (besselj(0,a)+b*bessely(0,a))*(1-T*a)-B*a*(besselj(1,a)+b*bessely(1,a))];
options = optimset('TolFun',1e-12,'TolX',1e-12,'Display','none');
alow = 0.01;
ahigh = 60.01;
blow = 0.01;
bhigh = 60.01;
da = 0.5;
db = 0.5;
a0 = alow:da:ahigh;
b0 = blow:db:bhigh;
icount = 0;
for i=1:numel(a0)
    for j=1:numel(b0)
        [x,r] = fsolve(@(x)fun(x(1),x(2)),[a0(i) b0(j)],options);
        if norm(r) < 1e-8
            icount = icount + 1;
            a(icount) = x(1);
            b(icount) = x(2);
            res(icount) = norm(r);
        end
    end
end
AB = [a(:),b(:),res(:)];
AB = sortrows(AB);
icount = 1;
M(1,:) = AB(1,:);
for i = 1:size(AB,1)-1
    if AB(i+1,1) - AB(i,1) > 0.1
        icount = icount + 1;
        M(icount,:) = AB(i+1,:);
    end
end
M

M = 12×3

0.0998 0.0039 0.0000 0.5303 1.6781 0.0000 6.6570 1.7032 0.0000 12.9507 1.7844 0.0000 19.2383 1.8170 0.0000 25.5238 1.8344 0.0000 31.8085 1.8452 0.0000 38.0926 1.8526 0.0000 44.3765 1.8579 0.0000 50.6603 1.8620 0.0000

Answer 2

[Torsten]

I don't know how the (an,bn) for your problem are enumerated, but maybe it's a start.

At least I can assure you that you will not be able to symbolically solve for an,bn the way you tried.

The equations are far too compliciated to allow analytical solutions.

T = 10;
B = 1;
K = 0.5;
fun = @(a,b)[(besselj(0,K*a)+b*bessely(0,K*a))*(1-T*a)+B*a*(besselj(1,K*a)+b*bessely(1,K*a));...
             (besselj(0,a)+b*bessely(0,a))*(1-T*a)-B*a*(besselj(1,a)+b*bessely(1,a))];
options = optimset('TolFun',1e-12,'TolX',1e-12)

options = struct with fields:
                   Display: []
               MaxFunEvals: []
                   MaxIter: []
                    TolFun: 1.0000e-12
                      TolX: 1.0000e-12
               FunValCheck: []
                 OutputFcn: []
                  PlotFcns: []
           ActiveConstrTol: []
                 Algorithm: []
    AlwaysHonorConstraints: []
           DerivativeCheck: []
               Diagnostics: []
             DiffMaxChange: []
             DiffMinChange: []
            FinDiffRelStep: []
               FinDiffType: []
         GoalsExactAchieve: []
                GradConstr: []
                   GradObj: []
                   HessFcn: []
                   Hessian: []
                  HessMult: []
               HessPattern: []
                HessUpdate: []
          InitBarrierParam: []
     InitTrustRegionRadius: []
                  Jacobian: []
                 JacobMult: []
              JacobPattern: []
                LargeScale: []
                  MaxNodes: []
                MaxPCGIter: []
             MaxProjCGIter: []
                MaxSQPIter: []
                   MaxTime: []
             MeritFunction: []
                 MinAbsMax: []
        NoStopIfFlatInfeas: []
            ObjectiveLimit: []
      PhaseOneTotalScaling: []
            Preconditioner: []
          PrecondBandWidth: []
            RelLineSrchBnd: []
    RelLineSrchBndDuration: []
              ScaleProblem: []
       SubproblemAlgorithm: []
                    TolCon: []
                 TolConSQP: []
                TolGradCon: []
                    TolPCG: []
                 TolProjCG: []
              TolProjCGAbs: []
                  TypicalX: []
               UseParallel: []

x = fsolve(@(x)fun(x(1),x(2)),[1 1],options);

Equation solved.

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

a = x(1)

a = 0.5303

b = x(2)

b = 1.6781

fun(a,b)

ans = 2×1

1.0e-15 * -0.8882 -0.4441

[Alex Sha]

The following results seem to be more precise：

a: 0.09975118859991
b: 0.00387367284443446

[mery]

@Alex Sha. How did you get this result?

a: 0.09975118859991?
b: 0.00387367284443446?

[Torsten]

How I can get the firs n-roots (an and bn), with n=20

This was exactly my questions: how do you want to enumerate/order your roots ?

According to a_n, starting from 0 ? And b_n follows somehow from your system ?

Or according to b_n, starting from 0 ? And a_n follows somehow from your system ?

Or another mechanism ?

[mery]

@Torsten, I need for example 20 first roots of a_n and 20 firs roots of b_n

(a_1,b_1); (a_2,b_2).............as a list,

maybe like a function, because in the sum I need to replace a_n and b_n by these values for n=1,2,3..... (so if I can do like the same code of besselj(0,x), I will use it in my sum)

I don't know if I explained my idea well or not

[Torsten]

Don't you understand what I mean ?

Say you got roots as

(1 4)

(2356,-34)

(-27.5,13.66)

(0.3+6*i,24-pi*i)

What is (a1,b1),(a2,b2),(a3,b3),(a4,b4) ?

You must define an ordering/enumeration somehow in order to insert them at the correct places in your infinite sum equation !

[mery]

@Torsten I mean by the first roos like that resul,

I serch a_n and b_n in interval

 0 <= a <= 60 && 0 <= b <= 60
out:{{a->0.0909091,b->0.063493},{a->0.386002,b->1.52548}
    ,{a->6.45491,b->1.36345},{a->12.7503,b->1.42685}
    ,{a->19.0383,b->1.45155},{a->25.324,b->1.46458}
    ,{a->31.6087,b->1.47262},{a->37.893,b->1.47807}
    ,{a->44.1769,b->1.482}
    ,{a->50.4607,b->1.48497},{a->56.7443,b->1.4873}} 

but I have no idea how I can enumerate them and insert them in my sum. That's why I need to define a function for a_n and b_n as x_k in besselj(0,x).

[Torsten]

And why does

a = 0.5303

b = 1.6781

not appear in the list ?

Did you choose different values for T, B and K in Mathematica ?

So in principle, the (an,bn) are ordered according to he size of an starting from 0 ?

[mery]

No its same B ,T and K

B = 1;
T = 10;
K = 0.5;

fun = {(BesselJ[0, K*a] + b*BesselY[0, K*a])*(1 - T*a) + 
   B*a*(BesselJ[1, K*a] + b*BesselY[1, K*a]), (BesselJ[0, a] + 
      b*BesselY[0, a])*(1 - T*a) - 
   B*a*(BesselJ[0, a] + b*BesselY[0, a])}

%FindRoot[fun, {a, 20}, {b, 20}]
NSolve[fun[[1]] == 0 && fun[[2]] == 0 && 0 <= a <= 60 && 0 <= b <= 60]

[Torsten]

Your function is incorrect (at least according to your equations in your original post).

fun = {(BesselJ[0, K*a] + b*BesselY[0, K*a])*(1 - T*a) - 
        B*a*(BesselJ[1, K*a] + b*BesselY[1, K*a]), (BesselJ[0, a] + 
        b*BesselY[0, a])*(1 - T*a) - 
        B*a*(BesselJ[1, a] + b*BesselY[1, a])}

instead of

fun = {(BesselJ[0, K*a] + b*BesselY[0, K*a])*(1 - T*a) - 
        B*a*(BesselJ[1, K*a] + b*BesselY[1, K*a]), (BesselJ[0, a] + 
        b*BesselY[0, a])*(1 - T*a) - 
        B*a*(BesselJ[0, a] + b*BesselY[0, a])}

[Torsten]

First you will have to find out in which order the (an,bn) are enumerated.

Then you have to find the first - say - 30 pairs in order and save them to a file. Since they won't change for the calculation of the infinite sum, you can do this step independent of the sum calculation.

Then you can program the sum calculation after reading in the pairs (an,bn) from file.

Since the step of finding the zeros can be done beforehand, I suggest you let MATHEMATICA calculate them and you just read them in your MATLAB code from file when needed.

[mery]

Yes, you are right, my original function is

fun = {(BesselJ[0, K*a] + b*BesselY[0, K*a])*(1 - T*a) + 
        B*a*(BesselJ[1, K*a] + b*BesselY[1, K*a]), (BesselJ[0, a] + 
        b*BesselY[0, a])*(1 - T*a) - 
        B*a*(BesselJ[0, a] + b*BesselY[0, a])}

I have an error in Mathematica

[Torsten]

So the system in your first post is incorrect ?

It looked nicely symmetric in construction ...

[mery]

fun = @(a,b)[(besselj(0,K*a)+b*bessely(0,K*a))*(1-T*a)+B*a*(besselj(1,K*a)+b*bessely(1,K*a));...

(besselj(0,a)+b*bessely(0,a))*(1-T*a)-B*a*(besselj(1,a)+b*bessely(1,a))];

this is my function like latex systeme, I have an error in Mathematica, but I still have to write this function

function ak,bk=findakbk(K, B, T, n,sstep)

[Torsten]

As I said: Copy the roots from MATHEMATICA. This is the easiest way.

[mery]

@Torsten thank you very much