Need help with this code- Chevyshev rational function on matlab
Ältere Kommentare anzeigen
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 Kommentare
John D'Errico
am 23 Dez. 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
am 23 Dez. 2014
John D'Errico
am 24 Dez. 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.
Antworten (0)
Kategorien
Mehr zu Polynomials finden Sie in Hilfe-Center und File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!