Método Compuesto de Simpson de Euler y Runge-Kutta de orden cuatro (Algoritmo y Código en Matlab)

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)

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.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&#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);

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

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