Simpsons, Euler y Runge-Kutta de orden cuatro

Esta es una vista previa del archivo. Inicie sesión para ver el archivo original

src/euler/euler.mfunction y=euler(n,a,b,h)
format long
%****************************************************************
%** Universidad Jose Antonio Paez **
%** Evaluación Final **
%** **
%** Alumna: Romina Betancourt ** 
%** Materia: Calculo Numerico **
%****************************************************************
fprintf('Universidad Jose Antonio Paez\n\Facultad de Ingenieria\n');
fprintf('Escuela de Ingenieria en Computacion\n');
fprintf('Calculo Numerico\n\n');
fprintf('Evaluacion Final\n\n');
fprintf('Alumna:\n');
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
src/euler/feuler.mfunction f=feuler(x,y)
f=-y+x+2;%es la funcion del ejemplo 1
src/runge-kutta/runge_kutta.mfunction runge_kutta
%****************************************************************
%** Universidad Jose Antonio Paez **
%** Evaluación Final **
%** **
%** Alumna: Romina Betancourt ** 
 %** Materia: Calculo Numerico **
%****************************************************************
clear;
clc;
%****************************************
fprintf('Universidad Jose Antonio Paez\n\Facultad de Ingenieria\n');
fprintf('Escuela de Ingenieria en Computacion\n');
fprintf('Calculo Numerico\n\n');
fprintf('Evaluacion Final\n\n');
fprintf('Alumna:\n');
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);
src/simpson/simpsons1.m%****************************************************************
%** Universidad Jose Antonio Paez **
%** Evaluación Final **
%** **
%** Alumna: Romina Betancourt ** 
 %** Materia: Calculo Numerico **
%****************************************************************
clear;
clc;
fprintf('Universidad Jose Antonio Paez\n\Facultad de Ingenieria\n');
fprintf('Escuela de Ingenieria en Computacion\n');
fprintf('Calculo Numerico\n\n');
fprintf('Evaluacion Final\n\n');
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)
doc/informe final.docxREPÚBLICA BOLIVARIANA DE VENEZUELA
UNIVERSIDAD JOSÉ ANTONIO PAEZ
ESCUELA DE INGENIERÍA EN COMPUTACIÓN
 (
Evaluación 
Práctica Final
)
Sección: 305C1
Romina Betancourt CI: 16052570
Catedra: Cálculo Numérico
San Diego, febrero, 2015
INDICE
PP
1. Método Compuesto de Simpson (6 puntos)
1.1) Algoritmo……………………………………………………………….….03
1.2) Programa.......................................................................................................03
2. Método de Euler
2.1) Algoritmo......................................................................................................05
2.2) Programa.......................................................................................................05
3. Método de Runge-Kutta de orden cuatro
3.1) Algoritmo......................................................................................................06
3.2) Programa.......................................................................................................06
Evaluación Final
1) Método Compuesto de Simpson (6 puntos)
1.3) 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)
XI1 = 0; (suma de f (x2i-1)
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)
Si no tomar XI1 = XI1 + f(x)
Paso 6: Tomar XI = h (XI0 + 2XI2 + 4XI1)/3
Paso 7: Salida (XI)
Parar
1.4) Programa: 
%****************************************************************
%** Universidad Jose Antonio Paez ** 
%** Evaluación Final ** 
%** **
%** Alumna: Romina Betancourt ** 
 %** Materia: Calculo Numerico **
%****************************************************************
clear;
clc;
fprintf('Universidad Jose Antonio Paez\n\Facultad de Ingenieria\n');
fprintf('Escuela de Ingenieria en Computacion\n');
fprintf('Calculo Numerico\n\n');
fprintf('Evaluacion Final\n\n');
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 (6 puntos)
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
%****************************************************************
%** Universidad Jose Antonio Paez **
%** Evaluación Final **
%** **
%** Alumna: Romina Betancourt ** 
%** Materia: Calculo Numerico **
%************************************************************
fprintf('Universidad Jose Antonio Paez\n\Facultad de Ingenieria\n');
fprintf('Escuela de Ingenieria en Computacion\n');
fprintf('Calculo Numerico\n\n');
fprintf('Evaluacion Final\n\n');
fprintf('Alumna:\n');
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 (8 puntos)
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
%****************************************************************
%** Universidad Jose Antonio Paez **
%** Evaluación Final **
%** **
%** Alumna: Romina Betancourt ** 
%** Materia: Calculo Numerico **
%****************************************************************
clear;
clc;
fprintf('Universidad Jose Antonio Paez\n\Facultad de Ingenieria\n');
fprintf('Escuela de Ingenieria en Computacion\n');
fprintf('Calculo Numerico\n\n');
fprintf('Evaluacion Final\n\n');
fprintf('Alumna:\n');
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&#160;&#160;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);
2
Informe prueba Final / Cálculo Numérico