MÉTODO DE HEUN INTRODUCCION Aplica Aplicarr método métodos s numéri numéricos cos para para aprox aproxima imarr soluci solucione ones s de alguna algunas s ecuaciones diferenciales, viendo así la importancia de los métodos numéricos que radica en la aparición de ecuaciones diferenciales que no pued pueden en reso esolver lvers se por mét método odos trad tradic icio iona nale les s, y de ahí 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 deter determin minaci ación ón de dos derivada derivadas s en el interv intervalo alo una una en el punto punto inicial y otra en el !nal". #as dos derivadas se promedian después con la !nalidad de o$tener una mejor estimación de la pendiente en todo el intervalo. %ste procedimiento, conocido como método de &eun, se presenta en forma gr'!ca en la siguiente !gura
Primera iteración del método Euler modifcado (ecordar (ecordaremos emos que, en el método de %uler, la pendiente al inicio. yi′ = f ( xi , yi )
...( 1)
se utili)a para extrapolar linealmente a yi +1 yi +1
= y0 + f ( xi , yi ) h
..... ( 2 )
%n el método est'ndar de %uler de$ería parar aquí. *in em$argo, en el métod método o de &eun &eun la yi +1 calculada en la ecuación +" no es la
respuesta !nal, sino una predicción intermedia. or consiguiente, la distinguimos con un superíndice -. #a ecuaci ación +" se llama ecuación simplemen mente te ecuación predicto predictora ra o simple predicto predictorr. a una estimación estimación de yi +1 que permite el c'lculo de una estimación de la pendiente al !nal del intervalo/ yi′+1
=
(
)
f xi +1 , yi +1 .... ( 3 )
Así, Así, se com$ com$in inan an las las dos dos pend pendien iente tes s 0ecu 0ecuac acio ione nes s 1" 1" y 2"3 2"3 para para o$tener una pendiente promedio en el intervalo/
yî′ + yi′+1 f ( xi , yi ) + f ( xi +1 , yi +1 ) ′ y = = 2
2
%st' st' pend pendie ient nte e pro promedio edio se util utili) i)a a despu espué és para ara linealmente desde yi hasta yi +1 con el método de %uler/ yi +1
= yi +
f ( xi , yi )
+ f ( xi +1 , yi +1 ) 2
extr xtrapol apolar ar
h .. .... ( 4 )
A esta ecuación 4" se conoce como ecuación correctora o simplemente corrector %l método de &eun es el único método predictor5corrector predictor5corrector de un solo paso.
redictor
→
6orrector
→
yi +1
= y0 + f ( xi , yi ) h
yi +1
= yi +
f ( xi , yi )
+ f ( xi +1 , yi +1 ) 2
h
7$serve que de$ido a que en la ecuación 4" aparece yi +1 a am$os lados del signo igual, entonces se puede aplicar en una forma iterativa. %s decir, una estimación anterior se utili)ar' de manera repetida para proporcionar una estimación mejorada de yi +1
ERROR RE!TI"O PORCENTU! ε t
=
yi j+1 − yij+−11 j i +1
y
E#ERCICIO$ RE$UETO$ E#ERCICIO %
×100%
yi j+−11
= iteracion anterior
yi j+1
= actual del corrector
(esolver el pro$lema de valor inicial 896:7( AU;U*:7 A(AA A
dy dx
=
x2 + 4
en [0,0. 0,0.25 25]]
y
y(0) (0) = 1
6on n>?
6alculamos h/
h=
b−a n
=
redictor 6orrector
0.25 − 0 5
= 0.05
→ →
yi0+1 = yi + f ( xi , yi ) h yi +1
= yi +
f ( xi , yi )
+ f ( xi +1 , yi0+1 ) 2
h
0 ara ara la sigu siguie ient nte e solu soluci ción ón prim primer ero o calc calcul ular arem emos os el pred predic icto torr yi +1
y
seguidamente llevaremos llevaremos este resultado a la ecuación del corrector yi +1 que vendría a ser la solución de &eun
rimera interacción/ i = 0
= 0 ; x1 = 0.05 ; y0 = 1 ; yi0+1 = yi + f ( xi , yi ) h 02 + 4 0 y1 = 1 + ÷ 0.05 1 y10 = 1.2 f ( xi , yi ) + f ( xi +1 , yi0+1 ) yi +1 = yi + h
x0
2
02 + 4 0.052 + 4 [ 1 ] + [ 1.2 ]÷ y1 = 1 + ÷ 0.05 2 ÷÷ y1 = 1.1834 segunda interacción/ i = 1
h = 0.05
= 0.05 ; x2 = 0.10 ; y0 = 1.1834 yi0+1 = yi + f ( xi , yi ) h 0.052 + 4 0 y2 = 1.1834 + ÷ 0.05 1.1834 0 y2 = 1.3525 0 f ( xi , yi ) + f ( xi +1 , yi +1 ) yi +1 = yi + h
x1
;
h = 0.05
2
0.052 + 4 0.10 2 + 4 [ 1 ] + [ 1.3525 ]÷ y2 = 1.1834 + ÷ 0.05 2 ÷÷ y2 = 1.3421 tercera interacción/ i = 2
= 0.10 ; x1 = 0.15 ; y2 = 1.3421 yi0+1 = yi + f ( xi , yi ) h 0.102 + 4 0 y3 = 1.3421 + ÷ 0 .0 5 1.3421 y30 = 1.4925 f ( xi , yi ) + f ( xi +1 , yi0+1 ) yi +1 = yi + h
x2
;
2
0.102 + 4 0.152 + 4 [ 1.3421 ] + [ 1.4925 ]÷ y3 = 1.3421 + ÷ 0.05 2 ÷÷ y3 = 1.4842
cuarta interacción/ i = 3
h = 0.05
= 0.15 ; x4 = 0.20 ; y3 = 1.4842 0 yi +1 = yi + f ( xi , yi ) h 0.152 + 4 0 y4 = 1.4842 + ÷ 0.05 1.4842 0 y4 = 1.6197 f ( xi , yi ) + f ( xi +1 , yi0+1 ) yi +1 = yi + h
x3
;
h = 0.05
2
0.152 + 4 0.20 2 + 4 [ 1.4842 ] + [ 1.6197 ]÷ y4 = 1.4842 + ÷ 0.05 2 ÷÷ y4 = 1.6143 quinta interacción/ i = 4
= 0.20 ; x5 = 0.25 ; y4 = 1.6143 yi0+1 = yi + f ( xi , yi ) h 0.202 + 4 0 y5 = 1.6143 + ÷ 0.05 1.6143 y50 = 1.7394 0 f ( xi , yi ) + f ( xi +1 , yi +1 ) yi +1 = yi + h
x4
;
h = 0.05
2
0.202 + 4 0.252 + 4 [ 1.6143 ] + [ 1.7394 ]÷ y5 = 1.6143 + ÷ 0.05 2 ÷÷ y5 = 1.7353 i
@predictivo
@correctivo
-
1.+---
1.124
1
1.2?+?
1.24+1
+
1.411?
1.44+
2
1.B1CD
1.B142
4
1.D2C4
1.D2?2
E#ERCICIO & dy
= −0.5 y + 4e0.8 x
x = 0 hata x = 4 con un dx tamaEo tamaEo de paso paso igual igual a 1. #a condic condición ión inicial inicial es en x = 0 ! y = 2 896:7( AU;U*:7 A(AA A
6on 6on el méto método do de &eun integ integre re
redictor 6orrector
→ →
yi0+1 yi +1
= yi + f ( xi , yi ) h = yi +
f ( xi , yi )
+ f ( xi +1 , yi0+1 ) 2
h
Antes de resolver el pro$lema numéricamente, se utili)a el c'lculo para determinar la siguiente solución analítica/ y = y = y = y = y = y =
4 1.3 4 1.3 4 1.3 4 1.3 4 1.3 4 1.3
(e0.8 x
− e−0.5 x ) + 2e−0.5 x
(e0.8 ( 0 )
− e −0.5( 0) ) + 2e−0.5(0 ) = 2
(e0.8 (1)
− e−0.5(1) ) + 2 e−0.5(1) = 6.1946314
(e0.8 ( 2 )
− e −0.5( 2) ) + 2e−0.5( 2) = 14.8439219
(e0.8 ( 3)
− e−0.5(3) ) + 2e−0.5(3) = 33.6771718
(e0.8 ( 4 )
− e −0.5( 4) ) + 2e−0.5( 4) = 75.3389626
0 ara ara la sigu siguie ient nte e solu soluci ción ón prim primer ero o calc calcul ular arem emos os el pred predic icto torr yi +1
y
seguidamente llevaremos llevaremos este resultado a la ecuación del corrector yi +1 que vendría a ser la solución de &eun i = 0
= 0 ; x1 = 1 ; y0 = 2 ; h = 1 yi0+1 = yi + f ( xi , yi ) h 0 0,8 ( 0 ) − 0.5(2))1 y1 = 2 + (4e y10 = 5 f ( xi , yi ) + f ( xi +1 , yi0+1 ) yi +1 = yi + h
x0
2
[4e0,8(0) − 0.5(2)] + [4e0,8(1) − 0.5(5)] y1 = 2 + ÷1 2 y1 = 6.7010819 cuando i = 1
= 1 ; x2 = 2 ; y1 = 6.7010819 ; h = 1 yi0+1 = yi + f ( xi , yi ) h 7010819 + (4 e 0 ,8(1) − 0.5(6.7010819))1 y20 = 6, 70 0 y2 = 12.2527 f ( xi , yi ) + f ( xi +1 , yi0+1 ) yi +1 = yi + h
x1
2
[4e0,8(1) − 0.5(6.7010819)] + [4 e0,8(2 ) − 0.5(12.2527)] y2 = 6.7010819 + ÷1 2 y2 = 16.3197819 cuando i = 2
= 2 ; x3 = 3 ; y2 = 16.3197819 ; h = 1 yi0+1 = y1 + f ( xi , yi ) h y30 = 16.3197819 + ( 4e 0,8 ( 2 ) − 0.5(16.31 97819))1 y30 = 27.9720 f ( xi , yi ) + f ( xi +1 , yi0+1 ) yi +1 = yi + h
x2
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
= 3 ; x4 = 4 ; y3 = 37.1992489 ; h = 1 yi0+1 = yi + f ( xi , yi ) h y40 = 37.1992489 + ( 4e 0 ,8 (3) − 0.5(37.19 92489))1 y40 = 62.6923 f ( xi , yi ) + f ( xi +1 , yi0+1 ) yi +1 = yi + h
x3
2
[4e 0,8(3) − 0.5(37.19924 89)] + [4e 0,8( 4 ) − 0.5(62.6923)] y4 = 37.1992489 + ÷1 2 y4 = 83.3377674
x 1 + 2 4
y verdadero +.-----B.1C4B21 4 14.42C+ 1C 22.BDD1D 1 D?.22CB +B
y&eun +.-----B.D-1-1 C 1B.21CD 1C 2D.1CC+4 C 2.22DDB D4
FetFG" -.-.1 C.C4 1-.4B 1-.B+
E#ERCICIO ' (C)ri*tian *ala* ccac)ura+ Un tanque cilíndrico de fondo plano con un di'metro de 1.?m contiene un = 1.5 Kg " L a una altura de 2m.*e desea sa$er la líquid líquido o de densid densidad ad ρ = altura del líquido dentro del tanque tres minutos después de que se a$re completamente la v'lvula de salida, la cual da un gasto de 0.6 A 2 ga #3 " s donde A es el 'rea seccional del tu$o de salida y es 78.5 ×10 −4 # 2 ! $9.81 #" 2 con B iteraciones.
$OUCION
G = 0.6 A 2 ga
Acumulación>entrada/
−0.6 A V
dV ρ dt
= 0
*alida/
2 ga
π
π
4
4
= ( 1.5) 2 a
( 1.5)
2
da dt
= −0.6 A
2 ga
e donde da dt
=−
2.4 A 2 ga π ( 1.5)
2
= −0.0026653
2 ga al considerar como tiempo cero el
momento de a$rir la v'lvula y adem's la altura $uscada a un tiempo de 1da
= −0.0026653
2 ga
dt s, se llega a a 0 ( )
= 3m a ( 180 ) = '
a%! t%&
tiliare#o *or el #etodo analitico
da
0.0026 2665 653 3 = −0.00
dt da
−
a
1 2
= 0.0026653 1
u% a
C
du =
2
=2
3
2 ga util utili iar are# e#o o e*a e*ara raci cion on de +ari +aria ale le 2 g dt
1 da 2
a
⇒ −∫
da a
= 0.0026653
∫
2 g dt
1
⇒ −2 a 2 + C = 0.0026653
nue nuetr traa ecua ecuaci cion on era era
al alla#o a#o t%18 t%180 0
1 2
2 a =
.44826 2682 823 3 ⇒ a ( 180) = 0.448
&allamos por el método de &eun
2 gt
- ( 0 ) = 0 .0026653 ⇒
3 − 0.0026653 2 gt g t 2
÷÷
2
2 g × 0 = 3
A. − yi +1 B. −
1 2
= yi + hf ( xi , yi ) → cao *redictor
× f ( xi , yi ) + f ( xi +1 , yi +1 )
C. − yi +1
1 = yi + h × f ( xi , yi ) + f ( xi +1 , yi +1 ) → cao corrector 2
i , % & ' . 1
t i
A
B.
C
2BC1+1?1-
+.2B?? 1.BDB 1.412D21.-+124 -.BC1? -.4+44+C
5-.-1C242 5-.-1D+4C 5-.-1?1?5-.-12-BB? 5-.-1-CD 5-.--D?B
+.41CD1 1.C?4 1.441C-1.-4DCD-.D1BD?
,/--0&10
/uetro error relati+o *orcentual era εt
=
yi j+1 − yij+−11 yi j+1
×100% ⇒ ε t =
0.44 0.4482 8268 6823 23 − 0.44 0.4482 8268 68 0.44826823
0 .0000513 ×100% = 0.
E#ERCICIO - (C)ri*tian *ala* ccac)ura+ Una solución de salmuera de salmuera ra)ón constante de B#Hmin, hacia el inte interio riorr de un depó depósi sito to que que inic inicial ialme ment nte e cont contie iene ne ?-# de solu soluci ción ón de salmuera en la cual se disolvieron GIg de sal. #a solución contenida en el depósito se mantiene $ien agitada y Juye hacia el exterior con la misma rapide). *i la concentración de sal presente en el depósito es de -.? IgH#, terminar la cantidad de sal presente en el depósito al ca$o de 1 minutos. K6uenta concentración alcan)ara de sal en el depósito en un tiempo de ? minL
$OUCION 6on h>1.+? y ? iteraciones
= kg de al dentro del de*oito en el intante t 3 dy y = 6 × 0.5 − 6 × = − 3 y ÷ 25 dt 50 y ( 0 ) = 5 y ( 5 ) = '
x ( t )
&allamos por el método de &eun A. − yi +1 B. −
1 2
= yi + hf ( xi , yi ) → cao *redictor
× f ( xi , yi ) + f ( xi +1 , yi +1 )
C. − yi +1
i , % & ' -
1 = yi + h × f ( xi , yi ) + f ( xi +1 , yi +1 ) → cao corrector 2
t i
A
B.
C
1.+? +.? 2.D? ?
1-.22D? 1+.2D1C 14.1+2C
+.++ 1.C14D? 1.B4C- 1.4+-+4
D.DD? 1-.1B? 1+.++22
%'/221%
DI!3R!M! DE FU#O
P$EUDOCODI3O
clear all dispMN%:77 % &%U=M" clc syms x syms y f>inlineinputMingrese f>inlineinputMingrese la derivada/M,MsM""O x>inputMingrese x>inputMingrese el valor de x/M"O y>inputMingrese el valor de y/M"O h>inputMingrese h>inputMingrese el valor de h/M"O n>inputMingrese n>inputMingrese numero de iteraciones/M"O clc dispMxi , yi , yiP1" , @iP1" M"O GyiP1"caso predictor G@iP1"caso corrector for i>1/n s>hPxO y1>fevalf,x,y"O hy1>hQy1O y+>yPhy1O y2>fevalf,s,y+"O hy+>y2QhO yn>yPhy1Phy+"H+"O fprintfMRnG-.1f G-.4f G-.4f G-.4fM,x ,y ,y+ ,yn "O y>ynO x>xPhO end
P!RTE &
E#ERCICIO$ PROPUE$TO$
1.5%n un tanque perfectamente agitado se tiene 4--# de una salmuera en la cual est'n disueltos de sal común =a6#", en cierto momento se hace llegar al tanque un gasto de un - salmuera que contiene -.? Sg de sal común por litro. *i tiene un gasto de salida de -#Hmin determine cri*tian *ala* ccac)ura " a" K
Re*pue*ta a)
dy
= 40 − 0.2 y
dx h1
! ( 0)
= 25
y ( 10 )
='
r*ta ! ( 10 )
= 175.947
b)a olucion e otiene hata la cantidad de al en el tanue no ca#ie con el tie#*o ! ( 50 )
= 199.991
+.5* +.5*e e hace hacen n reacc eaccio iona narr isot isotér érmi mica came ment nte e +B+B- g de acet acetat ato o de etil etilo o CH 3COOC2 H 5 con 1D?g de hidróx hidróxido ido de sodio NaOH en solución acuosa ajustando el volumen total a ? litros" para dar acetato de sodio CH 3COONa y el alcohol etílico C2H 5OH de acue acuerd rdo o con con la sigu siguie ient nte e ecua ecuaci ción ón estequiometria. cri*tian *ala* ccac)ura "
CH 3COOC2 H 5 + NaOH
→ CH 3COONa + C 2H 5OH
Re*pue*ta dy
0.59 − y ) ( 0.875 − y ) = 1.44 × 0.01 ( 0. dx h1
! (0)
=0
y ( 30 )
='
r*ta ! ( 30 )
= 0.169238
2.56al 2.56alcul cular ar el tiempo tiempo neces necesario ario para para que el nivel nivel del líquid líquido o dentr dentro o del tanque esférica con r>?m mostrando en la !gura pase de 4m a 2m.#a velocidad de salida por el ori!cio del fondo es de v>4.C?,el di'metro de dicho ori!cio es de 1-cm. cri*tian *ala* ccac)ura "
Re*pue*ta da dt
=−
0.012375 a
( 10a − a )
h10
2
a ( 0)
=4
a ( ')
=3
r*ta ! ( 1000 )
= 2.9796
Utili)aremos ahora el guide de ecuaciones diferenciales para poder desarrollar los ejercicios Utili)ando los ejercicicos propuestos
Pro4lema % #a parte a
Ahora desarrollaremos la parte $
Pro4lema &
Pro4lema '
desarro rolla llarr los siguie siguiente ntes s ejerc ejercici icios os 896 896:7 :7( ( AU AU;U ;U*: *:7 7 A(A A(AA A -/5desar A
dy dx
= x 2 − y
con y (1) = 2 ; calc calcule ule y (1.5) 1.5) ; *ara *ara h
= x + 2 y − 1
= 0.1
con y (2) (2) = 1 ; calcule y (2. (2.5) ; *ara h = 0.1
Corr Corre eremo remo* * el e6er 6ercici cicio o prop ropue* ue*to a+ en M!T M!T!7 896:7( AU;U*:7 A(AA A
%l resultado en y?>1.?C+
Corr Correr erem emo* o* el e6er e6erci cici cio o prop propue ue*t *to o 4+ en M!T M!T!7 !7 896:7( AU;U*:7 A(AA A
