Solve system for Roots of Bessel Function
9 views (last 30 days)
Show older comments
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
on 9 Aug 2022
In the infinite sum, only the an appear.
Could you give a literature link to this sum representation ?
Accepted Answer
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
0 Comments
More Answers (1)
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)
x = fsolve(@(x)fun(x(1),x(2)),[1 1],options);
a = x(1)
b = x(2)
fun(a,b)
15 Comments
See Also
Products
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!