%derivación 5 puntos
syms x
disp(\'Fórmula de 5 puntos\');
f=input(\'Escriba la función: \');
h=input(\'Escriba el valor de h: \');
x0=input(\'Escriba el valor de x0: \');
j=diff(f);
fprintf(\'Seleccione 1 o 2 según la fórmula que desea usar:\\n\\t1- Cinco términos \\n\\t2- Cuatro términos (fórmula centrada)\')
op=input(\'\\nopción: \');
if op==1
fprintf(\'\\n\')
x1=x0+h;
fprintf(\'x1= x0+h= %9.15f\\n\',double(x1))
x2=x0+2*h;
fprintf(\'x2= x0+2*h= %9.15f\\n\',double(x2))
x3=x0+3*h;
fprintf(\'x3= x0+3*h= %9.15f\\n\',double(x3))
x4=x0+4*h;
fprintf(\'x4= x0+4*h= %9.15f\\n\',double(x4))
yaprox=1/(12*h)*(-25*subs(f,x0)+48*subs(f,x1)-36*subs(f,x2)+16*subs(f,x3)-3*subs(f,x4));
fprintf(\'Valor aproximado: (1/12h)*(-25fx0+48fx1-36fx2+16fx3-3fx4)) = %9.15f\\n\',double(yaprox))
yex=subs(j,x0);
e=abs(yaprox-yex);
end
if op==2
fprintf(\'\\n\')
x1=x0-2*h;
fprintf(\'x1= x0-2*h= %9.15f\\n\',double(x1))
x2=x0-h;
fprintf(\'x2= x0-h= %9.15f\\n\',double(x2))
x3=x0+h;
fprintf(\'x3= x0+h= %9.15f\\n\',double(x3))
x4=x0+2*h;
fprintf(\'x4= x0+2*h= %9.15f\\n\',double(x4))
yaprox=1/(12*h)*(subs(f,x1)-8*subs(f,x2)+8*subs(f,x3)-subs(f,x4));
fprintf(\'Valor aproximado = (1/12h)(fx1-8fx2+8fx3-fx4)= %9.15f\\n\',double(yaprox))
yex=subs(j,x0);
e=abs(yaprox-yex);
end
fprintf(\'El valor exacto es: %9.15f\\n\',double(yex));
fprintf(\'El error es: %9.15f\\n\',double(e));
%Derivación 3 puntos________________________________________________________________________________________________________________
%derivación // Hecho por Diego Quintanilla
syms x
disp(\'Fórmula de 3 puntos\');
f=input(\'Escriba la función: \');
h=input(\'Escriba el valor de h: \');
x0=input(\'Escriba el valor de x0: \');
j=diff(f);
fprintf(\'Seleccione 1 o 2 según la fórmula que desea usar:\\n\\t1-Tres puntos con tres términos \\n\\t2-Tres puntos con dos términos (fórmula centrada)\')
op=input(\'\\nopción: \');
if op==1
fprintf(\'\\n\')
x1=x0+h;
fprintf(\'x1= x0+h = %9.15f\\n\',double(x1));
x2=x0+2*h;
fprintf(\'x2= x0+2*h = %9.15f\\n\',double(x2));
yaprox=1/(2*h)*(-3*subs(f,x0)+4*subs(f,x1)-subs(f,x2));
fprintf(\'Valor aproximado= (1/2h)(-3fx0+4fx1-fx2) = %9.15f\\n\',double(yaprox))
yex=subs(j,x0);
e=abs(yaprox-yex);
end
if op==2
fprintf(\'\\n\')
x1=x0+h;
fprintf(\'x1=x0+h = %9.15f\\n\',double(x1));
x2=x0-h;
fprintf(\'x1=x0-h = %9.15f\\n\',double(x2));
yaprox=1/(2*h)*(subs(f,x1)-subs(f,x2));
fprintf(\'El valor aproximado es: (1/2h)(fx1-fx2) = %9.15f\\n\',double(yaprox));
yex=subs(j,x0);
e=abs(yaprox-yex);
end
fprintf(\'El valor exacto es: %9.15f\\n\',double(yex));
fprintf(\'El error es: %9.15f\\n\',double(e));
%________________________________________________________________________________________________________________________________
%funcion entre comillas, o punto y coma al final
disp(\'1-->Extrapolacion con funcion\');
disp(\'2-->Extrapolacion con tabla\');
select = input(\'seleccion : \');
arreglo = zeros();
switch(select)
case 1
syms x
f=input(\'digite la funcion f(x) : \' );
N=input(\'digite la cantidad de iteraciones N(j) : \');
h=input(\'digite h : \');
x0=input(\'digite el punto inicial x0 : \');
fprintf(\'\\n\\n\');
for i=1: N
div=h;
d=1;
if(i~=1)
div = (h/2^(i-1));
d=2^(i-1);
end
superior = x0+div;
inferior = x0-div;
arreglo(1,i) = (1/(2*div))*(subs(f,superior)-subs(f,inferior));
fprintf(\'N(%1.0f)(h/%1.0f)=(1/(2h/%1.0f))*(F(x0+h/%1.0f)-F(x0-h/%1.0f))=%3.15f\\n\\n\',1,d,d,d,d,arreglo(1,i));
end
for i=1:N
m=1;
fprintf(\'\\n---------------------------------------------------------------------------------\\n\\n\');
for j=1:(N-i)
k=2^(j);
d=(4^(i)-1);
arreglo((i+1) , j)=(arreglo(i,(j+1)))+((arreglo(i,(j+1))-arreglo(i,j))/d);
fprintf(\'\\n\\n N(%1.0f)(h/%1.0f)=N(%1.0f)(h/%1.0f)+((N(%1.0f)(h/%1.0f)-N(%1.0f)(h/%1.0f))/%1.0f =%3.15f\',(i+1),m,i,k,i,k,(j+1),m,d,arreglo((i+1),j));
m=m*2;
Vaprox = arreglo((i+1),j);
end
end
syms x;
dif=diff(f,x);
Vreal = double(subs(dif , x0));
fprintf(\'\\nValor aproximado es =%9.15f \',Vaprox);
fprintf(\'\\nValor exacto es =%9.15f \' , Vreal );
err=abs(Vreal-Vaprox);
fprintf(\'\\nError =%3.3e \\n\\n\',err);
case 2
otherwise
return;
end
%n(2,2)=n(2,1)+[n(2,1)-n(1,1)]/k
%n(3,2)=n(3,1)+[n(3,1)-n(2,1)]/k
%n(4,2)=n(4,1)+[n(4,1)-n(3,1)]/k
%____________________________________________________________________________________________________________________________________
%Integracion de Romberg
format long
syms x y z;
g=input(\'Escriba el integrando g(x)=\');
n=input(\'Escriba el tamaño R(n,n) n: \');
a=input(\'limite a=\');
b=input(\'limite b=\');
RM=zeros(n,n);
fa=subs(g,a);
fb=subs(g,b);
for j=1:n
h(j)=(b-a)/(2^(j-1));
fprintf(\'Valores para h cuando j=1:n\')
h(j)
end
RM(1,1)=(h(1)/2)*(fa+fb);
for j=2:n
Sum=0;
for i=1:(2^(j-2))
c=(a+((2*i-1)*h(j)));
fc=subs(g,c);
Sum=Sum+fc;
end
RM(j,1)=(1/2)*(RM(j-1,1)+h(j-1)*Sum);
end
for j=2:n
for i=j:n
RM(i,j)=RM(i,j-1) + (RM(i,j-1)-RM(i-1,j-1))/((4^(j-1))-1);
end
end
for j=1:n
for i=j:n
if j==1&&i==1
fprintf(\'R(1,1)=h1/2[fa+fb]=%9.15f\\n\',RM(1,1))
elseif j==1&&i>1
fprintf(\'R(%1.0f,1)=(1/2)[R(%1.0f,1)+h%1.0f SumF(a+(2i-1)h%1.0f]=%9.15f\\n\',i,i-1,i-1,i,RM(i,1))
else
fprintf(\'R(%1.0f,%1.0f)=R(%1.0f,%1.0f)+R(%1.0f,%1.0f)-R(%1.0f,%1.0f)/(4^(%1.0f-1)-1)=%9.15f\\n\',i,j,i,j-1,i,j-1,i-1,j-1,j,RM(i,j))
end
end
end
fxd=input(\'Escriba la misma función pero en términos de y, f(y): \');
iex=int(fxd,y,a,b);
e=abs(RM(n,n)-iex);
fprintf(\'\\nValor aproximado: %9.15f\\n\',RM(n,n))
fprintf(\'\\nValor real: %9.15f\\n\',double(iex))
fprintf(\'\\nEl error es: %9.15f\\n\',double(e))
%_________________________________________________________________________________________________________________________
function D=nca(a,b,n) %Newton Cotes abierta
disp(\'Newton-Cotes abiertas\')
syms x
F=input(\'Escriba la función f(x): \');
f=inline(F);
h=(b-a)/(n+2);
fprintf(\'Seleccione la fórmula a aplicar: \\n\\t 1-Regla del punto medio (n=0) \\n\\t 2-Para n=1 \\n\\t 3-Para n=2 \\n\\t 4-Para n=3\\n\');
op=input(\'opción: \');
if op==1
x0=a+h;
fx0=f(x0);
inte=2*h*fx0;
iex=int(F,x,a,b);
e=abs(iex-inte);
end
if op==2
x0=a+h;
x1=x0+h;
fx0=f(x0);
fx1=f(x1);
inte=3*h/2*(fx0+fx1);
iex=int(F,x,a,b);
e=abs(iex-inte);
end
if op==3
x0=a+h;
x1=x0+h;
x2=x0+2*h;
fx0=f(x0);
fx1=f(x1);
fx2=f(x2);
inte=4*h/3*(2*fx0-fx1+2*fx2);
iex=int(F,x,a,b);
e=abs(iex-inte);
end
if op==4
x0=a+h;
x1=x0+h;
x2=x0+2*h;
x3=x0+3*h;
fx0=f(x0);
fx1=f(x1);
fx2=f(x2);
fx3=f(x3);
inte=5*h/24*(11*fx0+fx1+fx2+11*fx3);
iex=int(F,x,a,b);
e=abs(iex-inte);
end
fprintf(\'\\nEl valor aproximado es: %9.15f\',double(inte));
fprintf(\'\\nEl valor exacto es: %9.15f\',double(iex));
fprintf(\'\\nEl error es: %9.15f\\n\',double(e));
end
%______________________________________________________________________________________________________________________
function D=ncc(a,b,n)% Newton Cotes cerradas
disp(\'Newton-Cotes cerradas\')
syms x
F=input(\'Escriba la función f(x): \');
f=inline(F);
h=(b-a)/n;
fprintf(\'Seleccione la fórmula a aplicar: \\n\\t 1-Regla del trapecio \\n\\t 2-Regla de Simpson \\n\\t 3-Regla de 3/8 de Simpson \\n\\t 4-Regla de Boole\\n\');
op=input(\'opción: \');
if op==1
fa=f(a);
fb=f(b);
inte=h/2*(fa+fb);
iex=int(F,x,a,b);
e=abs(iex-inte);
end
if op==2
fa=f(a);
fb=f(b);
x1=a+h;
fx1=f(x1);
inte=h/3*(fa+4*fx1+fb);
iex=int(F,x,a,b);
e=abs(iex-inte);
end
if op==3
fa=f(a);
fb=f(b);
x1=a+h;
x2=a+2*h;
fx1=f(x1);
fx2=f(x2);
inte=3*h/8*(fa+3*fx1+3*fx2+fb);
iex=int(F,x,a,b);
e=abs(iex-inte);
end
if op==4
fa=f(a);
fb=f(b);
x1=a+h;
x2=a+2*h;
x3=a+3*h;
fx1=f(x1);
fx2=f(x2);
fx3=f(x3);
inte=2*h/45*(7*fa+32*fx1+12*fx2+32*fx3+7*fb);
iex=int(F,x,a,b);
e=abs(iex-inte);
end
fprintf(\'\\nEl valor aproximado es: %9.15f\',double(inte));
fprintf(\'\\nEl valor exacto es: %9.15f\',double(iex));
fprintf(\'\\nEl error es: %9.15f\\n\',double(e));
end
%____________________________________________________________________________________________________________________________________
%Compuesta de Simpson
fprintf(\'>>clc\\n\');
fprintf(\'>>syms x;\\n\');
n=input(\'>>n=\');
a=input(\'>>a=\');
b=input(\'>>b=\');
disp(\'>>h=(b-a)/n;\');
h=(b-a)/n;
disp(\'>>x0=a\');
a
for i=1:n
fprintf(\'>>x%1.0f=(a+%1.0f*h)\\n\',i,i);
z=a+(i)*h;
fprintf(\'\\n x%1.0f=\\n\\n\\t%1.15f\\n\\n\',i,z);
end
disp(\'>>x12=b\');
b
syms x;
g=input(\'>>g=\');
iex=int(g,x,a,b);
g
sum1=0;
sum2=0;
cad1=\'\';
cad2=\'\';
for i=1:n-1
y=a+i*h;
if(mod(i,2)==0)
sum1 = sum1 + subs(g,y);
nod=strcat(\'subs(g,x\' ,num2str(i) ,\')\');
cad1 = strcat(cad1 , \'+\' , nod);
else
nod2=strcat(\'subs(g,x\' ,num2str(i) ,\')\');
sum2= sum2 + subs(g,y);
cad2 = strcat(cad2 , \'+\' , nod2);
end
end
fprintf(\'>>sum1=%s \\n\',cad1);
sum1
fprintf(\'\\n\');
fprintf(\'>>sum2=%s \\n\',cad2);
sum2
fprintf(\'\\n\');
fa=subs(g,a);
disp(\'>>fa=subs(g,a);\');
fb=subs(g,b);
disp(\'>>fb=subs(g,b);\');
fprintf(\'>>faprox=(h/3)*(fa + 2*sum1 + 4*sum2 + fb)\\n\');
faprox=(h/3)* (fa + 2*sum1 + 4*sum2 + fb);
e=abs(iex-faprox);
fprintf(\'El valor aproximado es: %9.15f\\n\',double(faprox));
fprintf(\'El valor exacto es: %9.15f\\n\',double(iex));
fprintf(\'El error es: %9.15f\\n\',double(e));
%___________________________________________________________________________________________________________
%Trapecio compuesto:
function D=trapeciocomp(a,b,n)
syms y
f=input(\'Escriba la función f(y): \');
fa=subs(f,a);
fb=subs(f,b);
h=(b-a)/n;
iex=int(f,y,a,b);
%función suma = fs
fs=0;
for i=1:n-1;
x(i)=a+h*i;
fs=fs+subs(f,x(i));
end
inte=0.5*h*(fa+2*fs+fb);
e=abs(inte-iex);
fprintf(\'El valor aproximado de la integral es: %9.15f\\n\',double(inte));
fprintf(\'El valor exacto de la integral es: %9.15f\\n\',double(iex));
fprintf(\'El error es: %9.15f\\n\',double(e));
end
%_________________________________________________________________________________________________________________________
%Cuadratura Gaussiana
format long;
P=input(\'>> Escriba Polinomio como vector [1 0 ...]=\');
N=input(\'>> Escriba cuánto vale N=: \');
disp(\'>>r=roots(P)\');
r=roots(P); display(r); c=zeros(N);
syms x;
for i=1:N
fprintf(\'>>x%1.0f=%9.15f \\n\', i, r(i));
end
for i=1:N
num=1;
denom=1;
cad=\'\';
fprintf(\'>>c%1.0f=\',i);
for j=1:N
if(i~=j)
num= num * (x-r(j));
denom= denom * (r(i)-r(j));
a=[\'x\',num2str(i)];
b=[\'x\',num2str(j)];
if(j~= N)
cad=strcat(cad ,strcat(\'(\',a, \'-\' ,b , \')\' , \'*\'));
fprintf(\'(x-x%1.0f)*\', j);
else
fprintf(\'(x-x%1.0f)\', j);
cad=strcat(cad ,strcat(\'(\',a, \'-\' ,b , \');\'));
end
end
end
fprintf(\'/(%s)\', cad);
fprintf(\'\\n>>int(c%1.0f, -1 , 1 , x )\\n\',i);
c(i)=int((num/denom),x,-1,1);
c(i)
end
syms x t
%Respuestas del polinomio
f=input(\'Escriba el integrando: \');
pol=input(\'Escriba el polinomio f(x): \');
a=input(\'a= \');
b=input(\'b= \');
T=double(solve(pol,x));
n=length(T);
format long
%Integracion de polinomios de Lagrange
for i=1:n
M=1;
for j=1:n
if(i==j)
else
M=M*((x-T(j))/(T(i)-T(j)));
end
end
C(i)=double(int(M,x,-1,1));
end
%Sustitucion de la funcion integrando a la nueva funcion
V=(1/2)*(((b-a)*t)+a+b);
dV=(b-a)/2;
ft=subs(f,V)*dV;
S=0;
%Uso de las C y de la nueva funcion evaluada en las raices del polinomio
for i=1:n
S=S+(C(i)*subs(ft,T(i)));
end
double(S)
%Valor exacto
Exacto=(int(f,x,a,b));
fprintf(\'El valor aproximado es %9.15f \\n\',double(S));
fprintf(\'El valor exacto es %9.15f \\n\',double(Exacto));
fprintf(\'El error %e \\n\',double(abs(S-Exacto)));
%___________________________________________________________________________________________________________________________________
%Adaptivo de cuadratura
function D=adapcua(a,b)
syms x
f=input(\'Escriba la función f(x): \');
h=(b-a)/2;
fa=subs(f,a);
fb=subs(f,b);
fx1=subs(f,a+h/2);
fx2=subs(f,a+h);
fx3=subs(f,a+3*h/2);
iaprox=(h/6)*(fa+4*fx1+fx2)+(h/6)*(fx2+4*fx3+fb);
fprintf(\'El valor aproximado de la integral es: %9.15f\\n\',double(iaprox))
end
%__________________________________________________________________________________________________________________________________
%Segunda derivada
disp(\'Seguda derivada\')
f=input(\'Escriba la función: \');
syms x
h=input(\'Escriba el valor de h: \');
x0=input(\'Escriba el valor de x0: \');
sd=1/h^2*(subs(f,x0-h)-2*subs(f,x0)+subs(f,x0+h));
y=diff(f);
yy=diff(y);
y2=subs(yy,x0);
fprintf(\'El valor aproximado de la segunda derivada es: %9.15f\\n\',double(sd))
fprintf(\'El valor exacto de la segunda derivada es: %9.15f\\n\',double(y2));
%____________________________________________________________________________________________________________________________
%Compuesta de Simpson para integrales impropias singularidad límite
%izquierdo, cuando es derecho, hacer cambio x=z;
%Ahorita está diseñado para n=16;
function D=compimpr(a,b,n,p,sl)
h=(b-a)/n;
H=double(h)
x0=a;
%Aquí agregar o quitar x según sea el tamaño de n
x1=x0+h;x2=x1+h;x3=x2+h;x4=x3+h;x5=x4+h;x6=x5+h;x7=x6+h;x8=x7+h;
x9=x8+h;x10=x9+h;x11=x10+h;x12=x11+h;x13=x12+h;x14=x13+h;x15=x14+h;x16=x15+h;
syms x;
g=input(\'Función g(x): \');
d1=diff(g);d2=diff(d1);d3=diff(d2);d4=diff(d3);
g1=subs(g,sl); g2=subs(d1,sl); g3=subs(d2,sl); g4=subs(d3,sl); g5=subs(d4,sl);
px=g1+g2*(x-(sl))+g3*(x-(sl))^2/2+g4*(x-(sl))^3/6+g5*(x-(sl))^4/24;
f=(g-px)/(x-(sl))^p;
%Aquí agregar los valores de x
sum1=subs(f,x2)+subs(f,x4)+subs(f,x6)+subs(f,x8)+subs(f,x10)+subs(f,x12)+subs(f,x14);%suma los pares - el último
sum2=subs(f,x1)+subs(f,x3)+subs(f,x5)+subs(f,x7)+subs(f,x9)+subs(f,x11)+subs(f,x13)+subs(f,x15);%Suma todo impar
sumatoria1=double(sum1)
sumatoria2=double(sum2)
w=px/(x-(sl))^p;
%Sí el límite con sigularidad es el a (izquierdo), no evaluar subs(f,a) dentro de i1.
%Si es el b (derecho), no evaluar subs(f,b);
i1=double(h/3*(2*sum1+4*sum2+subs(f,b)))
i2=double(int(w,x,a,b))
iaprox=double(i1+i2);
fprintf(\'El valor aproximado es: %9.15f\\n\',double(iaprox))
end
%____________________________________________________________________________________________________________________________
%Derivadas
%función derivada
%asin(u)=u\'/sqrt(1-u^2)
%acos(u)=-u\'/sqrt(1-u^2)
%atan(u)=u\'/(1+u^2)
%acot(u)-u\'/(1+u^2)
%asec(u)=u\'/(abs(u)*sqrt(u^2 -1)
%acsc(u)=-u\'/(abs(u)*sqrt(u^2 -1)
%a^u=a^u(ln(u))u\'
%sinx = cosx
%cosx = -sinx
%tanx = sec^2
%cotx = -csc^2
%secx = secx tanx
%cscx = -cscx cotx
%integración
%sqrt(a^2-u^2) sqrt(a^2+u^2) sqrt(u^2-a^2)
%u=a sin u=a tan u= a sec
du= a cos du=a sec^2 du= a tan sec
sqrt = a cos sqrt= a sec sqrt= a tan
%1/sqrt(a^2-u^2) = arcosin(u/a)
%1/(a^2 +u^2)=(1/a)(arcotan(u/a))
%1/(u*sqrt(u^2 - a^2) = 1/a(arcosec(u/a))
%longitud de arco__________________________________
%fórmula=sqrt(1+^[f\'(x)]^2)
%area de superficie de revolución_____________________
%formula= 2pi*f(x)*sqrt(1+^[f\'(x)]^2)
%trabajo________________________________
%formula= (densidad)(volumen pequeño diferencial)(distancia)
%nota: bombeo es llenar, vaciar es distancia disco a bomba (k-y) llena es Y
%fuerza________________________________________
%formula= densidad*h(y)L(y)
%Arandela________________________________
%formula= pi(R^2 - r^2)
%eje a la izquierda o abajo eje derecha o arriba
%f(x) - eje eje -f(x)
%Capa cilindrica_____________________________
%formula= 2pi*p(x)h(x) nota: p(x)=x; h(x)= funcion de tope rec,tipico
%eje a la izquierda o abajo eje derecha o arriba
%f(x) - eje eje -f(x)
%método de discos_____________________________--
%formula= pi[f(x)]^2