Solve system for Roots of Bessel Function

17 views (last 30 days)
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
  3 Comments
Torsten
Torsten on 9 Aug 2022
In the infinite sum, only the an appear.
Could you give a literature link to this sum representation ?

Sign in to comment.

Accepted Answer

Torsten
Torsten on 10 Aug 2022
Edited: Torsten on 10 Aug 2022
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

More Answers (1)

Torsten
Torsten on 10 Aug 2022
Edited: Torsten on 10 Aug 2022
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
  15 Comments
mery
mery on 10 Aug 2022
Edited: mery on 10 Aug 2022
@Torsten thank you very much

Sign in to comment.

Community Treasure Hunt

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

Start Hunting!