1) Método
Compuesto de Simpson
1.1)
Algoritmo:
Entrada:
puntos extremos a; b; entero positivo m
Salida: Aproximación XI de la integral I
Paso 1: Tomar h = (b - a)/2m
Paso 2: Tomar XI0 = f(a) + f (b)
Salida: Aproximación XI de la integral I
Paso 1: Tomar h = (b - a)/2m
Paso 2: Tomar XI0 = f(a) + f (b)
XI1
= 0; (suma de f (x2i-1)
XI2 = 0; (suma de f (x2i))
XI2 = 0; (suma de f (x2i))
Paso
3: Para i = 1;…; 2m – 1
Paso 4: tomar x = a + ih
Paso5: Si i es par entonces tomar XI2 = XI2 + f(x)
Paso 4: tomar x = a + ih
Paso5: Si i es par entonces tomar XI2 = XI2 + f(x)
Si no tomar XI1 = XI1 + f(x)
Paso
6: Tomar XI = h (XI0 + 2XI2 + 4XI1)/3
Paso 7: Salida (XI)
Parar
Paso 7: Salida (XI)
Parar
1.2)
Programa:
%****************************************************************
%**
Romina Betancourt **
%****************************************************************
clear;
clc;
fprintf('Alumna:\n');
fprintf('Romina
Betancourt\n\n');
%Ingreso
de datos
disp('Regla
de simpson')
disp('Integrales
definidas')
fx=input('digite
la funcion f(x) = ','s');
%load
limiteinferior.txt
%a
= load('limiteinferior.txt');
%load
limitesuperior.txt
%b
= load('limitesuperior.txt');
%load
tol.txt
%tol
= load('tol.txt');
a=input('limite
inferior = ');
b=input('limite
superior = ');
tol=input('tolerancia
= ');
%Condiciones
iniciales
err(1)=100;ns=2;exito=0;
%Calculo
de la integral
while
exito==0
h=(b-a)/ns;
x=a:h:b;
y=eval(fx);
if
(rem(ns,3)==0) %simpson 3/8
Iaprox(ns-1)=3*h/8*(y(1)+y(ns+1)+3*sum(y(2:3:ns-1))+3*sum(y(3:3:ns))+2*sum(y(4:3:ns-2)));
elseif
(rem(ns,2)==0) %simpson 1/3
Iaprox(ns-1)=h/3*(y(1)+y(ns+1)+4*sum(y(2:2:ns))+2*sum(y(3:2:ns-1)));
else %
combinacion 3/8 + 1/3
Iaprox(ns-1)=h/3*(y(1)+y(ns-2)+4*sum(y(2:2:ns-3))+2*sum(y(3:2:ns-4)))+3*h/8*(y(ns-2)+3*y(ns-1)+3*y(ns)+y(ns+1));
end
if
ns>2 % calculo del error
err(ns-1)=abs((Iaprox(ns-1)-Iaprox(ns-2))/Iaprox(ns-1))*100;
if
err(ns-1)<tol
exito=1;
break;
end
end
ns=ns+1;
end
%Presentacion
de resultados
n=2:ns;
fprintf('n');
disp([' segm' '
integ' ' error'])
disp([n'
Iaprox' err' ]);
fprintf('n
se alcanzo la solucion con % g segmentos n',ns);
fprintf('n
la integral aproximada es %g n',Iaprox(ns-1));
plot(x,y)
2) Método de Euler
2.1)
Algoritmo: Este
algoritmo calcula la soluci´on del problema de valor inicial (1.1) en puntos
equidistantes x1 = x0 + h ,x2 = x0 + 2h, x3 = x0 + 3h, ・ ・ ・ , xN = x0 + Nh, aquí f es tal que (1.1) tiene una solución ´única en [x0, xN].
Entrada: Valores iníciales x0, y0,
tama˜no de paso h y n´umero de pasos N
Para n=0,...,N-1, hacer
xn+1 = xn + h
yn+1 = yn + hf(xn, yn)
Salida xn+1, yn+1
Parar
2.2)
Programa:
function y=euler(n,a,b,h)
format long
%****************************************************************
%**
Romina Betancourt **
%****************************************************************
fprintf('Romina Betancourt\n\n');
fprintf('\n ///// EULER/////\n')
x=a:h:n*h;
y=zeros(n,1);
y(1)=b;
for k=1:n
f=feuler(x(k),y(k));
y(k+1)=y(k)+h*f;
end
3) Método de Runge-Kutta de orden cuatro
3.1) Algoritmo: Este
algoritmo calcula la soluci´on del problema de valor inicial (1.1) en puntos
equidistantes x1 = x0 + h ,x2 = x0 + 2h, x3 = x0 + 3h, ・ ・ ・ , xN = x0 + Nh, aqu´ı f es tal
que (1.1) tiene una solución ´única en [x0, xN].
Entrada: Valores iníciales x0, y0,
tamaño de paso h y numero de pasos N
Para n=0,...,N-1, hacer
k1 = hf(xn, yn)
k2 = hf(xn + 1
2h, yn + 1
2k1)
k3 = hf(xn + 1
2h, yn + 1
2k2)
k4 = hf(xn + h, yn + k3)
yn+1 = yn + 1
6 (k1 + 2k2 + 2k3 + k4)
Salida xn+1, yn+1
Parar
3.2) Programa:
function
runge_kutta
%****************************************************************
%**
Romina Betancourt **
%****************************************************************
clear;
clc;
fprintf('Romina
Betancourt\n\n');
format
short % exportando formato
%
ingreso de Datos
fprintf('\n
///// RUNGE-KUTTA DE ORDEN 4/////\n')
dfx=input('\EC.
DIFERENCIAL dfx =','s');
a
=input('\n ingrese a=:\n');
b =input('\n
ingrese b=:\n');
%load limiteinferior.txt
%a = load('limiteinferior.txt');
%load
limitesuperior.txt
%b
= load('limitesuperior.txt');
y0=input('\n  EL
valos inicial y0 = :\n');
n=input('\n
cual es el No de Pasos? :\n');
h=(b-a)/n;
%vector
final agustado al numero de pasos
xs=a:h:b;
fprintf('\n''Funcion
descrita');
for
i=1:n
it=i-1;
a=xs(i);
x=a;
y =
y0;
k1=h*eval(dfx);
x=a+h/2;
y=y0+k1/2;
k2=h*eval(dfx);
x=a+h/2;
y=y0+k2/2;
k3=h*eval(dfx);
x=a+h;
y=y0+k3;
k4=h*eval(dfx);
y0=y0+(k1+2*k2+2*k3+k4)/6;
%formato
de impresion de resultados
fprintf('\n%2.0f%5.6f%10.6f\n',it,a,y0);
plot(xs,y0,':ok');
end
fprintf('\n
El punto aproximado y(x1) es = %8.6f\n',y0);
Archivo limiteinferior.txt (caso de leer desde archivo debe modificar el código, los valores adentro pueden ser los que usted guste)
-10
-9
-8
-7
-6
-5
-4
-3
-2
-1
0
1
2
3
4
5
6
7
8
9
10
-9
-8
-7
-6
-5
-4
-3
-2
-1
0
1
2
3
4
5
6
7
8
9
10
Archivo limitesuperior.txt (caso de leer desde archivo debe modificar el código, los valores adentro pueden ser los que usted guste)
10
9
8
7
6
5
4
3
2
1
0
-1
-2
-3
-4
-5
-6
-7
-8
-9
-10
9
8
7
6
5
4
3
2
1
0
-1
-2
-3
-4
-5
-6
-7
-8
-9
-10
