Michael Carey wrote:
> I have looked at linear regression for approximating a linear graph,
> but how would I approximate a square or third order polynomial?
>
> Thanks Mick C
>
real life example. Using the Peter's Eggleton stellar evolution code
I've got a table of the stellar interior parameters for the Brown
Dwarf (0.2 M_sun, 0.006 L_sun). The parameters are pressure,
temperature, luminosity, density and etc.. . Now I want to use their
functional form in my dynamo model. The best way is to get the
Chebyshev approximations to them. I do as follows,
1) read the table to maxima
2) find the spline interpolation
3) find the Chebyshev approximations
The code is below. The gz-table is in attachment.
best regards,
Valery
load("interpol");
load("numericalio.lisp");
load(orthopoly);
orthopoly_returns_intervals:false;
chebftGL(x1,x2,n,func):= block([xm,xl,fk,ck,numer],
xm:0.5*(x2-x1), xp:0.5*(x2+x1),fun:func,
xk:makelist(xp-xm*bfloat(cos(%pi*i/(n-1))), i, 0,n-1),
zk:makelist((xk[i+1]-xp)/xm, i,0,n-1),
fk:makelist(bfloat(fun(xk[i])),i,1,n),
ck[1]:float(sum(float(fk[i]),i,1,n)/(n-1)),
for j: 2 while j <= n do (
ck[j]:2./(n-1)*(bfloat(sum(bfloat(fk[i]*chebyshev_t(j-1,zk[i])),i,1,n))-
.5*(fk[1]*chebyshev_t(j-1,zk[1])+fk[n]*chebyshev_t(j-1,zk[n])))),
c:makelist(ck[i],i,1,n))$
model:read_matrix("fc.02.mdl1")$
G0:6.6726e-8;
MS:1.99e33;
RS:6.96e10;
g0:G0*MS/RS^2;
c_p:3.5e+08;
FS:6.34101764118e10;
r_pr:submatrix(model,2,4,5,6,7,8,9,10,11)$
r_m:submatrix(model,3,4,5,6,7,8,9,10,11)$
r_rh:submatrix(model,2,3,5,6,7,8,9,10,11)$
r_T:submatrix(model,2,3,4,6,7,8,9,10,11)$
r_l:submatrix(model,2,3,4,5,6,7,8,9,10)$
/* Firstly, we find the spline approximations */
define(press(x),cspline(r_pr))$
define(mass(x),cspline(r_m))$
define(rh(x),cspline(r_rh))$
define(temp(x),cspline(r_T))$
define(lum(x),cspline(r_l))$
xit:flatten(makelist(submatrix(r_pr,2)[i],i,1,length(r_pr)))$
xb:xit[3];
xe:xit[length(r_pr)-1];
chnod1(x):=(x-.5*(xe+xb))/(.5*(xe-xb))$
/* pressure vs R */
pr:chebftGL(xb,xe,13,press);
define(pre(x),sum(chebyshev_t(k,chnod1(x))*pr[k+1],k,0,12)-
5.007751721341311e+15)$
plot2d([[discrete,
transpose(r_pr)[1],transpose(r_pr)[2]],pre(x)],[x,xb,xe]);
/* density vs R */
rho:chebftGL(xb,xe,13,rh);
define(rhO(x),sum(chebyshev_t(k,chnod1(x))*rho[k+1],k,0,12)-5.645400818803552)$
plot2d([[discrete,
transpose(r_rh)[1],transpose(r_rh)[2]],rhO(x)],[x,xb,xe]);
/* temperature vs R */
tem1:chebftGL(xb,xe,13,temp);
define(tem(x),sum(chebyshev_t(k,chnod1(x))*tem1[k+1],k,0,12)-258851.0607949151)$
plot2d([[discrete,
transpose(r_T)[1],transpose(r_T)[2]],tem(x)],[x,xb,xe]);
/* mass vs R */
mas:chebftGL(xb,xe,13,mass);
define(massa(x),sum(chebyshev_t(k,chnod1(x))*mas[k+1],k,0,12)-.008322575912386927)$
plot2d([[discrete,
transpose(r_m)[1],transpose(r_m)[2]],massa(x)],[x,xb,xe]);
/* luminosity (in L_sun) vs R */
lum1:chebftGL(xb,xe,13,lum);
define(lumi(x),sum(chebyshev_t(k,chnod1(x))*lum1[k+1],k,0,12)-2.4089768666352515e-4)$
plot2d([[discrete,
transpose(r_l)[1],transpose(r_l)[2]],lumi(x)],[x,xb,xe]);
/* pressure stratification scale */
define(hp(x),(-pre(x)/diff(pre(x),x)))$
/* density stratification scale */
define(hrho(x),(rhO(x)/diff(rhO(x),x)))$
/* convective turn-over time [day] */
define(tau(x),(4*RS^4*c_p/FS/G0/MS*tem(x)*rhO(x)*hp(x)^2*x^4/lumi(x)/massa(x)/3)^(1/3))$
plot2d(tau(x)/3600/24,[x,xb,xe]);
> _______________________________________________
> Maxima mailing list
> Maxima at math.utexas.edu
> http://www.math.utexas.edu/mailman/listinfo/maxima
>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: fc.02.mdl1.gz
Type: application/x-gzip
Size: 6229 bytes
Desc: not available
Url : http://www.math.utexas.edu/pipermail/maxima/attachments/20081209/e5802bab/attachment.bin