%    integración de romber

clear all;

syms x y z q valorsitos;

disp('Metodo de Romberg')

f=input('Ingrese la funcion: ');

a=input('Ingrese el valor de a: ');

b=input('Ingrese el valor de b: ');

n=input('Ingrese el valor de N: ');

cont1=1;

%encontrando los valores de h

while(cont1<=n)

    h(cont1)=(b-a)/(2^(cont1-1));

    fprintf('h%d=%.9f\n\n',cont1,h(cont1));

    cont1=cont1+1;

end

%generando la primera columna

t(1,1)=(h(1)/2)*(subs(f,a)+subs(f,b));

fprintf('R(1,1)=h1/2*(f(a)+f(b))==>');

disp(t(1,1));

cont1=2;

sumatoria=0;

while(cont1<=n)

    cont2=1;

    fprintf('R(%1.0f,%1.0f)=',cont1,1);

    valorsitos(cont1,cont2)=q+h(cont1-1);

    while(cont2<=(2^(cont1-2)))

       

        valorsitos(cont1,cont2+1)=(q+a+(2*cont2-1)*h(cont1));

        sumatoria=sumatoria+subs(f,(a+(2*cont2-1)*h(cont1)));

        cont2=cont2+1;

    end

    t(cont1,1)=(1/2)*(t(cont1-1,1)+h(cont1-1)*sumatoria);

    fprintf(')==>');

    disp(t(cont1,1));

    sumatoria=0;

    cont1=cont1+1;

end

%generando el resto de la tabla

cont1=2;

cont2=2;

cont3=2;

while(cont2<=n)

    while(cont1<=n)

       t(cont1,cont2)=t(cont1,cont2-1)+(t(cont1,cont2-1)-t(cont1-1,cont2-1))/(4^(cont2-1)-1);

       fprintf('R(%1.0f,%1.0f)=R(%1.0f,%1.0f)+(R(%1.0f,%1.0f)-R(%1.0f,%1.0f)/%d)==>',cont1,cont2,cont1,cont2-1,cont1,cont2-1,cont1-1,cont2-1,(4^(cont2-1)-1));

       disp(t(cont1,cont2));

       cont1=cont1+1;

    end

    cont3=cont3+1;

    cont1=cont3;

    cont2=cont2+1;

end

%clear all            integracion de simpsom

syms x y;

funcion=input('funcion=');

a=input('a=');

b=input('b=');

n=input('n=');

h=(b-a)/n;

for i=1:n-1

    x(i)=a+i*h;

    fprintf('x(%d)= ',i)

    pretty(x(i));

    %disp(x(i));

end

part2=0;

for i=1:(n/2-1)

    part2=part2+subs(funcion,x(i*2));

end

part3=0;

for i=1:(n/2)

    part3=part3+subs(funcion,x(i*2-1));

end

resultado=(h/3)*(subs(funcion,a)+2*part2+4*part3+subs(funcion,b));

double(resultado)

                           % extrapolación de richarson programa 1

function N = richardson2B(f, x0, n, h, pre)

%meter lo de arriba  asi como va copy paste! y  la funcion f es  primero

%hayq  declararla!

% f es la funcion, x0 el de aproximar, n el subindice de N, h es h y pre es

% precisi�n (opcional)

if nargin < 5

    pre = 8;

end

pre = num2str(pre);

N = zeros(n,n);

% generamos la primera columna

for a = 1:n

    h_ = h/(2^(a-1));   

    expr = sprintf('[1/(2*(%g))]*[f(%g) - f(%g)]', h_, x0 + h_, x0 - h_);  

    N(a, 1) = centrada3_fun(f, x0, h_, str2num(pre));   

   

end

% las demas columnas

for a = 2:n

    for b = a:n

        N(b, a) = N(b, a-1) + (N(b, a-1) - N(b-1, a-1))/(4^(a-1) - 1);

    end

end

                      % extrpolacion de Richardson programa 2

function RICHARDSON(h,N,Xo,g)

%h es el espaciado, N es hasta q numero se va a llegar N3,N4, etc. Xo es el

%punto en donde se va a evaluar la derivada y g es la funcion.

R=zeros(N,N);

for j=1:N

    X1=Xo+h;

    X2=Xo-h;

    syms x

    g1=subs(g,X1);

    g2=subs(g,X2);

    %Primero usamos la formula de 3 puntos

    %f'(h)=(1/2h)(f(x0+h)-f(x0-h))

    %f'(h/2)=(1/2(h/2))(f(x0+(h/2))-f(x0-(h/2)))

    R(j,1)=(1/(2*h))*(g1-g2);

    fprintf('h=%1.4f\n',h)

    h=h/2;

end

for j=2:N

    for i=j:N

        R(i,j)=R(i,j-1) + (R(i,j-1)-R(i-1,j-1))/((4^(j-1))-1);

    end

end

for j=1:N

   for i=j:N

fprintf('N(%d,%d)= %3.10f\n',i,j,R(i,j))

   end

end

valor_aproximado=R(N,N)

gp=diff(g);

valor_exacto=subs(gp,Xo)

error=abs(valor_aproximado - valor_exacto)

%Formula de Richardson

%               N1(h/2)-N1(h)       

%N2(h)=N1(h/2)+ -------------

%                4^(j-1))-1

%

%                 N1(h/4)-N1(h/2)       

%N2(h/2)=N1(h/4)+ -------------

%                  4^(j-1))-1

%

%....