Need help with this code- Chevyshev rational function on matlab
Show older comments
I am trying to approximate some data. Matlab couldn't approximate polynomial so I approximated using another software which can generate matlab file. This is the actual function.
%F(x,y)=(a+dT1(x )+eT1(y )+hT2(x)+iT2(y)+lT3(x )+mT3(y )+pT4(x)+qT4(y )+tT5(x )+uT5(y )... +abT6(x)+acT6(y)+afT7(x)+agT7(y )+ajT8(x )+akT8(y )+anT9(x)+aoT9(y)+arT10(x)+... asT10(y))/((1+bT1(x)+cT1(y)+fT2(x)+gT2(y )+jT3(x)+kT3(y )+nT4(x)+oT4(y )+... rT5(x )+sT5(y)+vT6(x)+aaT6(y)+adT7(x)+aeT7(y)+ahT8(x)+aiT8(y)+alT9(x)+amT9(y 1)+... apT10(x)+aqT10(y ) )+atT11(x )+auT11(y ))
% where Tn(x) = 2xTn-1(x) - Tn-2(x), n≥2 is chebysheb rational function
I do understand in this code is I need to provide input x and y and matrix c is the coefficient of above function. I get lost in function z = evalcratl(order, logx, logy, x, y, p,... s0, s1, s2, s3, s4, s5, s6, s7, s8, s9). How is this function working. Meaning of Evalcratl and tcnt .
Any help on this would be appreciated.
--------------------------------------------------------------------------------
[rowx colx]=size(xa);
if(rowx~=1 & colx~=1)
error('x must be scalar or 1D array');
return;
end
[rowy coly]=size(ya);
if(rowy~=1 & coly~=1)
error('y must be scalar or 1D array');
return;
end
c=[
0.7823326556648382,
-0.003073888389053987,
6.311956873270135E-16,
-0.0003235701036352635,
0,
0.1184529306672700,
0.7000618410087535,
-0.06699883001305451,
0.6882427553087074,
-0.001289845952369463,
7.834476942891433E-16,
0.0003263428678678636,
0,
-0.01433571833230648,
-0.2173515431820155,
0.01552171410016730,
-0.2211421352264687,
-0.0001289154504598550,
6.559509737139252E-16,
0.0004529420876541368,
0,
-0.01218726621626016,
-0.09562258140958195,
0.003063157214620869,
-0.08732570079084124,
0.0002432826238712180,
4.144106621110242E-16,
0.0003501722099832917,
0,
-0.008018886525623193,
-0.04998847489715176,
0.003499439465963643,
-0.04330683933474295,
6.716713718394032E-05,
1.191551227290313E-16,
8.914863235537785E-05,
0,
0.001270750287501435,
-0.0001942757598005615,
];
lenx=length(xa);
leny=length(ya);
for(j=1:leny)
for(i=1:lenx)
x=xa(i);
y=ya(j);
z=evalcratl(38,0,0,x,y,c,...
0.000000000000000,5.000000000000000,...
0.000000000000000,0.000000000000000,...
0.000000000000000,5.000000000000000,...
0.000000000000000,0.000000000000000,...
0.1962500000000000,0.8042500000000000);
za(i,j)=z;
end
end
%--------------------------------------------------------------
function z = evalcratl(order, logx, logy, x, y, p,...
s0, s1, s2, s3, s4, s5, s6, s7, s8, s9)
%--------------------------------------------------------------
tx=[];
ty=[];
if(logx==0)
x=(x-s0)/s1;
else
x=(log(x)-s2)/s3;
end
if(logy==0)
y=(y-s4)/s5;
else
y=(log(y)-s6)/s7;
end
if (order==6)
tcnt=3;
elseif (order==10)
tcnt=4;
elseif (order==14)
tcnt=5;
elseif (order==18)
tcnt=6;
elseif (order==22)
tcnt=7;
elseif (order==26)
tcnt=8;
elseif (order==30)
tcnt=9;
elseif (order==34)
tcnt=10;
elseif (order==38)
tcnt=11;
elseif (order==42)
tcnt=12;
else
return;
end
if(tcnt>7)
if(x<-1.0)
x=-1.0;
end
if(x>1.0)
x= 1.0;
end
if(y<-1.0)
y=-1.0;
end
if(y>1.0)
y= 1.0;
end
end
tx(1)=1.0;
ty(1)=1.0;
tx(2)=x;
ty(2)=y;
for j=2:1:tcnt-1
tx(j+1)=2*x*tx(j)-tx(j-1);
ty(j+1)=2*y*ty(j)-ty(j-1);
end
m=2;
num=p(1);
den=1.0+p(2)*tx(m)+p(3)*ty(m);
for(j=3:4:order-1)
num=num+p(j+1)*tx(m);
num=num+p(j+2)*ty(m);
m=m+1;
den=den+p(j+3)*tx(m);
den=den+p(j+4)*ty(m);
end
if (den==0)
z=0;
else
z=(num/den)*s8+s9;
end
return;
Thank You in advance
3 Comments
John D'Errico
on 23 Dec 2014
I'm pretty sure that matlab COULD have solved for that approximating polynomial, just that you did not know how to do it. There is a big difference.
sagar pokhrel
on 23 Dec 2014
John D'Errico
on 24 Dec 2014
Again, those limits were more about your understanding of the numerical analysis involved than about the true capabilities of a numerical tool.
As far as an explanation of the code you have been given, I don't see much purpose to it. Sorry, but reading undocumented code that does something you don't understand is a losing proposition. Either trust the provider or don't use it.
Answers (0)
Categories
Find more on Polynomials in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!