PoliOrtogonales2.html

> restart;with(numapprox):Digits :=3;

Se puede expandir una funci—n polin—mica en tŽrmino de polinomios ortogonales

> poli := 1+2*x+3*x^3 +4*x^5;
a := chebyshev(poli,x);

> restart;with(OrthogonalSeries) :

Igualmente, ese mismo polinomio puede ser expresado en varias bases de polinomios ortogonales

> poli := 1+2*x+3*x^3 +4*x^5;

> S1 := ChangeBasis(poli,ChebyshevT(n,x)) ;

> S2 := ChangeBasis(poli,HermiteH(n,x)) ;

> S3 := ChangeBasis(S1,HermiteH(n,x));

Los polinomios ortogonales son funciones y como tales se evalœan

> LegendreP(2,x);
LegendreP(2,0.3);simplify(%);

igualmente se expresan como series hipergeomŽtricas

> LegendreP(a,z):
% = convert( %, hypergeom);

> ?convert;

> restart;f(x):=exp(-(x -1/2)^2);

Si suponemos una funci—n arbitraria se calculan los coeficientes

> a[n]:=int(LegendreP(n,x)*f(x),x=-1..1)/int(LegendreP(n,x)^2,x=-1..1);

> N := 5;
for n from 0 to N do
a[n]:=evalf( int(LegendreP(n,x)*f(x),x=-1..1) ) / evalf( int(LegendreP(n,x)^2,x=-1..1)) ;
end do;

> S(x):= sum(a[j]*LegendreP(j,x),j=0..N);

> ST(x):= taylor( f(x), x=0, 7+1 );
STP(x):= convert(ST(x), polynom):

> plot([[x,S(x),x=-1..1],[x,STP(x),x=-1..1],[x,f(x),x=-1..1]],color=[red,blue,black]);

Interpolando puntos a funciones polin—micas...

> restart;with(plots):
XYdatos:= [[2, 8], [4,10], [6,11], [8,18], [10,20], [12,34]]: # o equivalentemente con solo puntos
listplot(XYdatos,style = POINT, symbol = DIAMOND, color=red);

Warning, the name changecoords has been redefined

se reescalan al intervalo ('-1,1)

> Npuntos := 6;

> for i from 1 to Npuntos do
XYdatosRE[i,1]:= XYdatos[i,1]/5 +(5-12)/5:
XYdatosRE[i,2]:= XYdatos[i,2]:
end do:

Se plantea un sistema de ecuaciones a partir de la expansi—n del vector como combinaci—n lineal de los seis primeros polinomios de Legendre

> for i from 1 to Npuntos do
EcuacP[i] := simplify(-XYdatosRE[i,2]+ sum(C[j]*LegendreP(j,XYdatosRE[i,1]),j=0..(Npuntos-1))=0 );
end do;

Se resuelve el sistema para encontrar las componentes del vector en el subespacio 6 dimensional

> solve({EcuacP[1],EcuacP[2],EcuacP[3],EcuacP[4],EcuacP[5],EcuacP[6]},{C[0],C[1],C[2],C[3],C[4],C[5]} );assign(%);

> for i from 1 to Npuntos do
EcuacT[i] := simplify(-XYdatosRE[i,2]+ sum(A[j]*ChebyshevT(j,XYdatosRE[i,1]),j=0..(Npuntos-1))=0 );
end do;

> solve({EcuacT[1],EcuacT[2],EcuacT[3],EcuacT[4],EcuacT[5],EcuacT[6]},{A[0],A[1],A[2],A[3],A[4],A[5]} );assign(%);

con lo cual el polinomio de interpolaci—n ser‡

> PInterpolP(x):= evalf(sum(C[j]*LegendreP(j,x),j=0..(Npuntos-1)));

> PInterpolT(x):= evalf(sum(A[j]*ChebyshevT(j,x),j=0..(Npuntos-1)));

> GrafInterpolP:= plot(PInterpolP(x),x=-1..1):
GrafInterpolT:= plot(PInterpolT(x),x=-1..1):

> GrafPuntos:=listplot([[XYdatosRE[1,1],XYdatosRE[1,2]],
                [XYdatosRE[2,1],XYdatosRE[2,2]],
                [XYdatosRE[3,1],XYdatosRE[3,2]],
                [XYdatosRE[4,1],XYdatosRE[4,2]],
                [XYdatosRE[5,1],XYdatosRE[5,2]],
                [XYdatosRE[6,1],XYdatosRE[6,2]]],
                style = POINT, symbol = DIAMOND, color=blue):

se grafican los puntos y su polinomo de interpolacion

> display({GrafInterpolP,GrafPuntos}, title = "Puntos e Interpolacion Legendre");

> display({GrafInterpolT,GrafPuntos}, title = "Puntos e Interpolacion Tchebyshev");

Cuadratura de Gauss

Mostraremos la forma que toma la cuadratura de Gauss para n puntos. La cuadratura de Gauss permite expresar una integral como una sumatoria pesada del integrando evaluado en puntos estrat'egicos que son las ra'ices de los polinomios de Legendre.

> restart;

> N := 6;
PNx:=simplify(LegendreP(N,x));

> x := [fsolve(PNx)];

> Cuadratura := f -> sum(p['i']*f(x['i']),'i'=1..N);

> Ecuacs := []:

> for j from 1 to N do
Ecuacs := [op(Ecuacs),
Cuadratura(z -> z^(j-1)) = int(z^(j-1),z=-1..1)]:
end do:
Ecuacs;

> SistEc := {op(Ecuacs)}:
variab := {seq(p[k],k=1..N)};

> pesos := solve(SistEc,variab);assign(%);

> PoliGenerico2N1 := randpoly(z,dense,degree=2*N-1);

> fpoli := unapply(PoliGenerico2N1,z);

> IntegCuad:= Cuadratura(fpoli);

> IntegExact := int(PoliGenerico2N1,z=-1..1);evalf(%);

> errorCuad:= evalf(IntegExact - IntegCuad);

> PoliGenerico2N := randpoly(z,dense,degree=2*N);

> fpoli := unapply(PoliGenerico2N,z);

> IntegCuad:= Cuadratura(fpoli);

> IntegExact := int(PoliGenerico2N,z=-1..1);evalf(%);

> errorCuad:= evalf(IntegExact - IntegCuad);

> ff:=z->exp(-(z -1/2)^2);

> IntegCuad:= Cuadratura(ff);

> IntegExact := int(ff(z),z=-1..1);evalf(%);

> errorCuad:= evalf(IntegExact - IntegCuad);