Polynominal approximation, least squares fit & Bessel function

iamfromspacebaby

Hi all,

As Bessel functions are not supported in Max/MSP I have to do a polynominal approximation...

I'm using this webpage: //oak.ucc.nau.edu/jws8/dpgraph/drumheads.html as a resource.

The approximation is p(r) = rm (1 - r2) (a + b r2 + c r4 ) where 0 < r < 1. It is an approximation for Jm(km,n r).

Apparently I have to do a "least squares fit" for a, b and c.

Where m = 2, n = 1 and K = 5.135 the author of the webpage gets a = 3.03 and b = -2.45. However the values I get are alot smaller, they are: a = -0.042855 and b = 0.065353.

The code I have written is below:

function [p]=calcPolyApproximation(m,K)
%n = order, K=m'th zero of the n'th order bessel function
bessel = besselj(m,(0:.01:1).*K);
[p,s] = polyfit((0:.01:1).*K,bessel,3);
z = polyval(p,(0:0.01:1).*K);
figure;plot(bessel,'*');
hold on;plot(z,'g');

Any help would be appreciated,

Carl.

Btw. I am very new to Matlab so sensitive with me

shwaip

expanding the polynomial:
p(r) = r^m (1 - r^2) (a + b*r^2 + c*r^4)
means you actually need a r^(m+2+4) term in your polyfit, and you're only using the first 3.

iamfromspacebaby

shwaip wrote:

expanding the polynomial:
p(r) = r^m (1 - r^2) (a + b*r^2 + c*r^4)
means you actually need a r^(m+2+4) term in your polyfit, and you're only using the first 3.

Ahh yes, but even if I change that I still get very small results? Any other ideas?

shwaip

you can simplify this a bunch if you're willing to assume that 'm' is even. If this is ok, you can change the polynomial you're fitting. Let me know if this is an acceptable assumption - I'll whip up some code quickly.

shwaip

actually, you don't need that assumption. I'll have something in a sec...

shwaip

First you say:

p(r) = rm (1 - r2) (a + b r2 + c r4 ) where 0 < r < 1.

but then you try to fit this bessel function:
bessel = besselj(m,(0:.01:1).*K);

which has r outside the range of 0,1.

I think this is probably part of the problem...

iamfromspacebaby

shwaip wrote:

First you say:

p(r) = rm (1 - r2) (a + b r2 + c r4 ) where 0 < r < 1.

but then you try to fit this bessel function:
bessel = besselj(m,(0:.01:1).*K);

which has r outside the range of 0,1.

I think this is probably part of the problem...

Thanks for the quick reply. Well I have included *K because in the formula it stated that K needs to be multipled by r... tbh I wasn't sure about that myself.

iamfromspacebaby

Thanks for your help so far... I think it is better but not 100%, and can someone enlighten me on the -01 -03 etc and how to get rid of them?

function [p]=calcPolyApproximation(m,K)
%n = order, K=m'th zero of the n'th order bessel function
bessel = besselj(m,0:.01:1,K);
[p,s] = polyfit(0:.01:1,bessel,4);
z = polyval(p,0:.01:1);
figure;plot(bessel,'*');
hold on;plot(z,'g');

iamfromspacebaby

Well, infact I dont think the above is using the K value.. damn.

If anyone else has more suggestions I'd appreciate it as I'm stumped.