jueves, 20 de octubre de 2011

newto divididas

function newton(x,b,xk)
n=length(xk);
for i=1:n
   f(i,1)=b(i);
end
for j=2:n
   k=j-1;
   for i=1:n+1-j
      f(i,j)=(f(i+1,k)-f(i,k))/(a(i+j-1)-a(i));
   end
end
prod=1;
sum=0;
for j=1:n
   sum=sum+f(1,j)*prod;
   prod=prod*(xk-x(j));
end

fprintf('\n        x           F[]\n');
fprintf('=============');
for j=1:n
   fprintf('=============');
end
fprintf('\n');
for i=1:n
   fprintf('%10.6f',x(i));
   for j=1:n+1-i
      fprintf('%10.6f',f(i,j));
   end
end
P=sum;
disp('EL RESULTADO DE LA INTERPOLACION ES:');
disp(P);

newton infinitas

function infinitas(x,y,xb)
m(:,1)=x;
m(:,2)=y;
c=0;
r=zeros(length(x)-1,length(x)+1);
for i=3:length(x)+1
    c=c+1;
    for j=1:length(y)-c
        m(j,i)=y(j+1)-y(j);
    end
   y=m(:,i);
end
p(1)=1;
p1(1)=1;
c=length(x);
for i=1:length(x)-1
    for j=1:i
        p1(2)=-x(j);
        p=conv(p,p1);
    end
    k=1;
    in=c;
    ind2=1;
    while k<=length(p);   
         r(i,in)=p(ind2);
            k=k+1;
            in=in+1;
            ind2=ind2+1;
    end
    c=c-1;
    p=ones(1);
end
h=x(2)-x(1);
disp(m);
ind=1;
for i=1:length(x)+1
    for j=3:length(x)+1
         r(ind,i)=r(ind,i)*(m(1,j)/(fact(ind)*(h^ind)));
         ind=ind+1;
    end
    ind=1;
end
pol=zeros(1,length(x)+1);
for i=1:length(x)+1
    for j=1:length(x)-1
        pol(i)=pol(i)+r(j,i);
    end
end
pol(length(x)+1)=pol(length(x)+1)+m(1,2);
f=poli(pol,xb);
disp(r);
disp(pol);
disp(f);
-------------------------------------------------------- 
function s=poli(x,y)
s=0;
n=length(x);
for i=1:n
s=s+x(i)*y^(n-i);
end
-----------------------------------------------------
function f=fact(x)
f=1;
for i=1:x
f=f*i;
end

miércoles, 12 de octubre de 2011

Lagrange

function Lagrange(x,y,xk)
n=length(x);
pol1(1)=1;
for i=1:n
    p1=1;
    p2=1;
    for j=1:n
        if(j~=i)
            p1=p1*(xk-x(j));
            p2=p2*(x(i)-x(j));
        end
    end
    d(i)=y(i)/p2;
    l(i)=p1/p2;
end
c(1)=1;
 for i=1:n
     z=length(c);
    
     for j=1:n
        if(j~=i)
            pol1(2)=-x(j);
            c=conv(c,pol1);
        end
     end
     for a=1:length(c)
         m(i,a)=c(a);
     end
     c=ones(z);
 end

 for i=1:length(m)
     for j=1:n
     m(j,i)=m(j,i)*d(j);
     end
 end
 for i=1:length(m)
     s=0;
     for j=1:n
         s=s+m(j,i);
     end
     r(i)=s;
 end
     p=0;
for i=1:n
    p=p+l(i)*y(i);
end
s=0;
n=length(r);
for i=1:n
    s=s+r(i)*(xk^(n-i));
end
 disp 'EL RESULTADO ES:';
disp(r);
disp(s);

EL RESULTADO ES:
     0     1     0     0
    6.7600

polinomio

function x=polinomio(a,b,xb)
n=length(a);
m=zeros(n);
p=n;
for i=1:n
    for j=1:n
        m(i,j)= (a(i)^(p-1));
        p=p-1;
    end
    p=n;
end
x=inv(m)*b';
p=n;
r=0;
for i=1:n
    r=r+x(i)*(xb^(p-1));
    p=p-1;
end
disp('los coef del polinomio')
disp(x)
disp('el resultado es de p(x)=')
disp(r);

por ejemplo si a=[0 1 2 3] y b=[0 1 4 9] y xb=2.6
el resul est

los coef del polinomio
   -0.0000
    1.0000
   -0.0000
         0
el resultado es de p(xb)=
    6.7600