Metodo De Heun(modificado) 38281z

MÉTODO DE HEUN INTRODUCCION Aplicar métodos numéricos para aproximar soluciones de algunas ecuaciones diferenciales, viendo así la importancia de los métodos numéricos que radica en la aparición de ecuaciones diferenciales que no pueden resolverse por métodos tradicionales, y de ahí la necesidad de implementar algún método de aproximación.

DEFINICION Un método para mejorar la estimación de la pendiente emplea la determinación de dos derivadas en el intervalo (una en el punto inicial y otra en el final). Las dos derivadas se promedian después con la finalidad de obtener una mejor estimación de la pendiente en todo el intervalo. Este procedimiento, conocido como método de Heun, se presenta en forma gráfica en la siguiente figura

Primera iteración del método Euler modificado Recordaremos que, en el método de Euler, la pendiente al inicio.

yi  f  xi , yi 

...  1

se utiliza para extrapolar linealmente a yi 1

yi 1  y0  f  xi , yi  h .....  2  En el método estándar de Euler debería parar aquí. Sin embargo, en el método de Heun la yi 1 calculada en la ecuación (2) no es la

respuesta final, sino una predicción intermedia. Por consiguiente, la distinguimos con un superíndice 0. La ecuación (2) se llama ecuación predictora o simplemente predictor. Da una estimación de yi 1 que permite el cálculo de una estimación de la pendiente al final del intervalo:





yi1  f xi 1 , yi 1 ....  3

Así, se combinan las dos pendientes [ecuaciones (1) y (3)] para obtener una pendiente promedio en el intervalo:



yî  yi1 f  xi , yi   f xi 1 , yi 1  2 2

y 



Está pendiente promedio se utiliza después para linealmente desde yi hasta yi 1 con el método de Euler: yi 1  yi 



f  xi , yi   f xi 1 , yi 1 2

h

extrapolar

....  4 

A esta ecuación (4) se conoce como ecuación correctora o simplemente corrector El método de Heun es el único método predictor-corrector de un solo paso.

Predictor



yi 1  y0  f  xi , yi  h

Corrector



yi 1  yi 



f  xi , yi   f xi 1 , yi 1 2

h

Observe que debido a que en la ecuación (4) aparece yi 1 a ambos lados del signo igual, entonces se puede aplicar en una forma iterativa. Es decir, una estimación anterior se utilizará de manera repetida para proporcionar una estimación mejorada de yi 1 ERROR RELATIVO PORCENTUAL

t 

yij1  yij11 100% yij1

yij11  iteracion anterior yij1  actual del corrector

EJERCICIOS RESUELTOS EJERCICIO 1

Resolver el problema de valor inicial (VICTOR AUGUSTO ARAPA AQUINO)

dy x 2  4  dx y

en [0,0.25]

y (0)  1

Con n=5

Calculamos h:

h

b  a 0.25  0   0.05 n 5

Predictor



Corrector



yi01  yi  f  xi , yi  h yi 1  yi 

f  xi , yi   f  xi 1 , yi01  2

h

Para la siguiente solución primero calcularemos el predictor

yi01

y

seguidamente llevaremos este resultado a la ecuación del corrector yi 1 que vendría a ser la solución de Heun

Primera interacción:

i0 x0  0

;

x1  0.05

;

y0  1 ;

yi01  yi  f  xi , yi  h  02  4  0.05  1 

y10  1   y10  1.2

yi 1  yi 

f  xi , yi   f  xi 1 , yi01 

h 2  02  4 0.052  4  ][ ]  [ 1 1.2  0.05 y1  1   2       y1  1.1834

segunda interacción: i  1

h  0.05

x1  0.05 ;

x2  0.10

;

y0  1.1834

;

h  0.05

yi01  yi  f  xi , yi  h  0.052  4 y  1.1834    0.05  1.1834  0 2

y20  1.3525 yi 1  yi 

f  xi , yi   f  xi 1 , yi01   

y2  1.1834     

h 2 0.052  4 0.102  4  [ ][ ] 1 1.3525  0.05 2   

y2  1.3421 tercera interacción:

i2

x2  0.10

x1  0.15

;

;

y2  1.3421 ;

yi01  yi  f  xi , yi  h  0.102  4 y  1.3421    0.05  1.3421  0 3

y30  1.4925 yi 1  yi 

f  xi , yi   f  xi 1 , yi01   

y3  1.3421     

h 2 0.102  4 0.152  4  [ ] [ ] 1.3421 1.4925  0.05 2   

y3  1.4842

cuarta interacción: i  3

h  0.05

x3  0.15

x4  0.20

;

y3  1.4842

;

;

h  0.05

yi01  yi  f  xi , yi  h  0.152  4 y  1.4842    0.05  1.4842  0 4

y40  1.6197 yi 1  yi 

f  xi , yi   f  xi 1 , yi01  

h 2 0.152  4 0.202  4  [ ][ ] 1.4842 1.6197  0.05 2   



y4  1.4842      y4  1.6143

quinta interacción: i  4

x4  0.20

x5  0.25

;

; y4  1.6143 ; h  0.05

yi01  yi  f  xi , yi  h  0.202  4  0.05  1.6143 

y50  1.6143   y50  1.7394 yi 1  yi 

f  xi , yi   f  xi 1 , yi01   

y5  1.6143     

h 2 0.202  4 0.252  4  [ ][ ] 1.6143 1.7394  0.05 2   

y5  1.7353

i

Ypredictivo

Ycorrectivo

0

1.2000

1.1834

1

1.3525

1.3421

2

1.4115

1.4842

3

1.6197

1.6143

4

1.7394

1.7353

EJERCICIO 2

dy  0.5 y  4e0.8 x x  0 hasta x  4 con un dx tamaño de paso igual a 1. La condición inicial es en x  0 y y  2 (VICTOR Con el método de Heun integre AUGUSTO ARAPA AQUINO)

Predictor



Corrector



yi01  yi  f  xi , yi  h yi 1  yi 

f  xi , yi   f  xi 1 , yi01  2

h

Antes de resolver el problema numéricamente, se utiliza el cálculo para determinar la siguiente solución analítica:

4 0.8 x 0.5 x (e  e )  2e0.5 x 1.3 4 0.8(0) 0.5(0) y (e e )  2e0.5(0)  2 1.3 4 0.8(1) 0.5(1) y (e e )  2e0.5(1)  6.1946314 1.3 4 0.8(2) 0.5(2) y (e e )  2e0.5(2)  14.8439219 1.3 4 0.8(3) 0.5(3) y (e e )  2e 0.5(3)  33.6771718 1.3 4 0.8(4) 0.5(4) y (e e )  2e0.5(4)  75.3389626 1.3 y

Para la siguiente solución primero calcularemos el predictor

yi01

y

seguidamente llevaremos este resultado a la ecuación del corrector yi 1 que vendría a ser la solución de Heun

i0 x0  0

x1  1 ; y0  2

;

;

h 1

yi01  yi  f  xi , yi  h y10  2  (4e0,8(0)  0.5(2))1 y10  5 yi 1  yi  

y1  2  

f  xi , yi   f  xi 1 , yi01 

h 2 [4e0,8(0)  0.5(2)]  [4e0,8(1)  0.5(5)] 1 2 

 y1  6.7010819 cuando i  1

x1  1 ;

x2  2

y1  6.7010819

;

h 1

;

yi01  yi  f  xi , yi  h y20  6, 7010819  (4e0,8(1)  0.5(6.7010819))1 y20  12.2527 yi 1  yi 

f  xi , yi   f  xi 1 , yi01 

h 2  [4e0,8(1)  0.5(6.7010819)]  [4e0,8(2)  0.5(12.2527)] y2  6.7010819   1 2   y2  16.3197819

cuando i  2

x2  2

;

x3  3

;

y2  16.3197819

;

h 1

yi01  y1  f  xi , yi  h y30  16.3197819  (4e0,8(2)  0.5(16.3197819))1 y30  27.9720 yi 1  yi 

f  xi , yi   f  xi 1 , yi01 

h 2  [4e0,8(2)  0.5(16.3197819)]  [4e 0,8(3)  0.5(27.9720)] y3  16.3197819   1 2   y3  37.1992489 cuando i  3

x3  3 0 i 1

y

;

x4  4

;

y3  37.1992489

;

h 1

 yi  f  xi , yi  h

y40  37.1992489  (4e0,8(3)  0.5(37.1992489))1 y40  62.6923 yi 1  yi 

f  xi , yi   f  xi 1 , yi01  

y4  37.1992489   

y4  83.3377674

h 2 [4e 0,8(3)  0.5(37.1992489)]  [4e0,8(4)  0.5(62.6923)] 1 2 

y verdadero 2.000000 0 6.194631 4 14.84392 19 33.67717 18 75.33896 26

x 0 1 2 3 4

yHeun 2.000000 0 6.701081 9 16.31978 19 37.19924 89 83.33776 74

|et|(%) 0.00 8.18 9.94 10.46 10.62

EJERCICIO 3 (Christian salas ccachura) Un tanque cilíndrico de fondo plano con un diámetro de 1.5m contiene un líquido de densidad   1.5 Kg / L a una altura de 3m.Se desea saber la altura del líquido dentro del tanque tres minutos después de que se abre completamente la válvula de salida, la cual da un gasto de 0.6 A 2 ga m3 / s donde

78.5 10

A 4

es 2

el

área 2

seccional

del

m y g=9.81 m/s con 6 iteraciones.

SOLUCION

tubo

de

salida

y

es

G  0.6 A 2 ga

Acumulación=entrada:

dV  0 dt

Salida:

0.6 A 2 ga

V

 2 da  0.6 A 2 ga  1.5 4 dt

 2  1.5  a 4

De donde

2.4 A 2 ga da   0.0026653 2 ga al considerar como tiempo cero el 2 dt   1.5  momento de abrir la válvula y además la altura buscada a un tiempo de 180

da  0.0026653 2 ga dt s, se llega a a 0  3m  

a=y t=x

a  180   ?

Utilizaremos por el metodo analitico da  0.0026653 2 ga utilizaremos separacion de variables dt da da  1  0.0026653 2 gdt    1  0.0026653 2 g  dt a2 a2 1 1 1 da u= a 2 du   2 a 2  C  0.0026653 2 gt  C  0   0 .0026653 2 g  0  3 2 a  2 3  0.0026653 2 gt C 2 3 nuestra ecuacion sera a      2   Hallamos t=180  a  180   0.44826823 Hallamos por el método de Heun

2

A.  yi 1  yi  hf  xi , yi   caso predictor





1 B.    f  xi , yi   f xi 1 , yi 1   2   1  C.  yi 1  yi  h    f  xi , yi   f xi 1 , yi 1      2 





 caso corrector

i

ti

A

B.

C

0 1 2 3 4 5 6

0 30 60 90 120 150 180

2.38655 1.86876 1.413730 1.02134 0.691588 0.424429

-0.019343 -0.017249 -0.015150 -0.0130665 -0.01097 -0.0088756

2.41971 1.89854 1.441900 1.047970 0.716758 0.448268

Nuestro error relativo porcentual sera yij1  yij11 0.44826823  0.448268 t  100%   t   100%  0.0000513 j yi 1 0.44826823

EJERCICIO 4 (Christian salas ccachura) Una solución de salmuera de salmuera razón constante de 6L/min, hacia el interior de un depósito que inicialmente contiene 50L de solución de salmuera en la cual se disolvieron %kg de sal. La solución contenida en el depósito se mantiene bien agitada y fluye hacia el exterior con la misma rapidez. Si la concentración de sal presente en el depósito es de 0.5 kg/L, terminar la cantidad de sal presente en el depósito al cabo de 1 minutos. ¿Cuenta concentración alcanzara de sal en el depósito en un tiempo de 5 min?

SOLUCION Con h=1.25 y 5 iteraciones

x  t   kg de sal dentro del deposito en el instante t dy 3  y  6  0.5  6     3  y dt 25  50 y  0  5 y  5  ?

Hallamos por el método de Heun

A.  yi 1  yi  hf  xi , yi   caso predictor





1 B.    f  xi , yi   f xi 1 , yi 1   2   1  C.  yi 1  yi  h    f  xi , yi   f xi 1 , yi 1      2 





 caso corrector

i

ti

A

B.

C

0 1 2 3 4

0 1.25 2.5 3.75 5

8 10.3375 12.37189 14.1239

2.22 1.91475 1.64908 1.420284

7.775 10.165 12.2233 13.9961

DIAGRAMA DE FLUJO

PSEUDOCODIGO

clear all disp('METODO DE HEUN') clc syms x syms y f=inline(input('ingrese la derivada:','s')); x=input('ingrese el valor de x:'); y=input('ingrese el valor de y:'); h=input('ingrese el valor de h:'); n=input('ingrese numero de iteraciones:'); clc disp('xi , yi , y(i+1) , Y(i+1) '); %y(i+1)caso predictor %Y(i+1)caso corrector for i=1:n s=h+x; y1=feval(f,x,y); hy1=h*y1; y2=y+hy1; y3=feval(f,s,y2); hy2=y3*h; yn=y+((hy1+hy2)/2); fprintf('\n%0.1f %0.4f %0.4f %0.4f',x ,y ,y2 ,yn ); y=yn; x=x+h; end

PARTE 2 EJERCICIOS PROPUESTOS 1.-En un tanque perfectamente agitado se tiene 400L de una salmuera en la cual están disueltos de sal común (NaCL), en cierto momento se hace llegar al tanque un gasto de un 80 salmuera que contiene 0.5 Kg de sal común por litro. Si tiene un gasto de salida de 80L/min determine (cristian salas ccachura) a) ¿Qué cantidad de sal hay en el tanque transcurrido 10 minutos? b) ¿Qué cantidad de sal hay en el tanque transcurrido un tiempo muy grande? Respuesta

a)

dy  40  0.2 y dx h=1

y  0   25 y  10   ?

rpta: y  10   175.947

b)La solucion se obtiene hasta la cantidad de sal en el tanque no cambie con el tiempo y  50   199.991 2.-Se hacen reaccionar isotérmicamente 260 g de acetato de etilo

CH 3COOC2 H 5 con 175g de hidróxido de sodio NaOH en solución acuosa (ajustando el volumen total a 5 litros) para dar acetato de sodio CH 3COONa y el alcohol etílico C2 H 5OH de acuerdo con la siguiente ecuación estequiometria. (cristian salas ccachura)

CH 3COOC2 H 5  NaOH  CH 3COONa  C2 H 5OH Respuesta

dy  1.44  0.01 0.59  y   0.875  y  dx h=1

y  0  0 y  30   ?

rpta: y  30   0.169238

3.-Calcular el tiempo necesario para que el nivel del líquido dentro del tanque esférica con r=5m mostrando en la figura pase de 4m a 3m.La velocidad de salida por el orificio del fondo es de v=4.895,el diámetro de dicho orificio es de 10cm. (cristian salas ccachura)

Respuesta

da 0.012375 a  dt  10a  a 2  h=10

a  0  4 a  ?  3

rpta: y  1000   2.9796

Utilizaremos ahora el guide de ecuaciones diferenciales para poder desarrollar los ejercicios Utilizando los ejercicicos propuestos Problema 1 La parte a

Ahora desarrollaremos la parte b

Problema 2

Problema 3

4.-desarrollar

los

siguientes

ejercicios

(VICTOR

AUGUSTO

ARAPA

AQUINO)

dy  x2  y dx

con y (1)  2 ; calcule y (1.5) ; para h  0.1

dy  x  2 y  1 con y (2)  1 ; calcule y(2.5) ; para h  0.1 dx

Correremos el ejercicio propuesto a) en MATLAB (VICTOR AUGUSTO ARAPA AQUINO)

El resultado en y5=1.8592 Correremos el ejercicio propuesto b) en MATLAB (VICTOR AUGUSTO ARAPA AQUINO)