TRABAJO DE INVESTIGACIÓN
MÉTODOS NUMÉRICOS EN ECUACIONES DIFERENCIALES ORDINARIAS
INTRODUCCIÓN
Las Las ecuaci ecuacion ones es dife difere renci ncial ales es sirv sirven en para para mode modela larr prol prolem emas as de cien cienci cias as e in!enier"a #ue re#uieren el camio de una variale respec$o a o$ra% En la ma&or par$e de ellos 'a& #ue resolver un prolema de valor inicial( es decir( resolver una ecuaci)n diferencial #ue sa$isface una condici)n inicial dada% En !eneral( las si$uacion si$uaciones es de la vida real( la ecuaci)n ecuaci)n diferencial diferencial #ue modela el prolema resul$a demasiado complicada para resolverla con una ma&or precisi)n( por lo #ue se recurre a dos $*cnicas para apro+imar la soluci)n% La primera consis$e en simplificar la ecuaci)n diferencial de modo #ue se pueda resolver precisa precisamen$ men$e e & u$ili, u$ili,ar ar despu* despu*ss la soluci soluci)n )n de la ecuaci) ecuaci)n n simpli simplific ficada ada para apro+imar la soluci)n de la ecuaci)n ori!inal% La se!unda( se vale de m*$odos para apro+imar la soluci)n del prolema ori!inal% Es$a $*cnica es la #ue se emplea por lo re!ular( pues los m*$odos de apro+imaci)n dan resul$ados m-s e+ac$os & una informaci)n realis$a sore el error% En la presen$e inves$i!aci)n se evidenciar- los dis$in$os m*$odos empleados para la resoluci)n de ecuaciones diferenciales ordinarias%
Trabajo Trabajo de Investigación Investigación ~Métodos Numéricos~
Página 1
Integración de Romberg Es$e m*$odo usa dos es$imaciones es$imaciones de una in$e!ral para calcular una $ercera m-s e+ac$a% b
Sea l ( h ) el valor de la in$e!ral #ue apro+ima a I =∫a f ( ( x ) dx ( median$e una
par$icipaci)n de suin$ervalos de lon!i$ud
h=
b −a ( usando la re!la de $rapecio% n
En$onces( el valor de la in$e!ral e+ac$a es la suma del valor es$imado con un paso h m-s el error !enerado. I = I ( h ) + E ( h )
Dond Donde e
E ( h ) es el error de $runcamien$o #ue se come$e al aplicar la re!la
$rapecial% El m*$odo de Romer! propone una $*cnica para acelerar la conver!encia( en el c-lculo c-lculo de una apro+imaci)n apro+imaci)n de una in$e!raci)n in$e!raci)n num*rica( para o$ener un $ercer valor m-s e+ac$o% La in$e!raci)n de Romer! se convier$e en el al!ori$mo m-s eficien$e den$ro de *s$e m*$odo( la cual es una f)rmula recursiva%
Trabajo Trabajo de Investigación Investigación ~Métodos Numéricos~
Página 2
Suponiendo #ue se $ienen dos apro+imaciones.
I ( h1)
e
I ( h2)
( con
suin$ervalos h1 & h2 respec$ivamen$e% I = I ( h1)+ E ( h1) I = I ( h2 ) + E ( h 2)
}
(h1 ) + E ( h1 )= I ( h2 ) + E (h 2)
⇒ I
( 1)
El error de la aplicaci)n m/l$iple de la re!la $rape,oidal se o$iene por. E ≅−
Donde
b− a 2 h f (ε) 12
f (ε)
es un promedio de la dole derivada en$re cier$os valores #ue
per$enecen a cada uno de los suin$ervalos% Si se despe0a E ( h 1) se puede conocer el error sin conocer f .
( )
h1 E ( h 1) ≅ E ( h2) h2
2
Sus$i$u&endo en la i!ualdad de las dos apro+imaciones( $enemos #ue.
()
2
I ( h 1) + E ( h2 )
h1 ≈ I (h2 ) + E ( h2 ) h2
I ( h 1) − I ( h2 ) ≈ E ( h2 )− E ( h 2)
() h1 h2
2
[ ( )]
h ¿ E ( h 2) 1− 1 h2
2
Trabajo de Investigación ~Métodos Numéricos~
Página 3
1or ende se lo!ra despe0ar E (h 2) . E ( h 2) ≈
I ( h1 )− I ( h2 )
( )
1−
h1
2
h2
I ( h1 )− I ( h2 ) I = I ( h 2) + 2 h1 1− h2
()
2a en el caso especial cuando
h2=
h1 2
siendo es$e el al!ori$mo de Romer!(
$enemos. I (h1 )− I (h2 ) I ≈ I ( h2 ) + 2 1−2 I ( h1) 4 I ≈ I ( h 2) − 3 3
1ara lo!rar en$ender el m*$odo( es convenien$e pensar #ue se $raa0a en niveles de apro+imaci)n% En un primer nivel( es cuando aplicamos la re!la del $rapecio( & para poder usar la f)rmula an$erior( deemos de duplicar cada ve, el n/mero de suin$ervalos. as" podemos comen,ar con un suin$ervalo( lue!o con dos( cua$ro( oc'o( e$c%( 'as$a donde se pida% Lue!o( se pasa al se!undo nivel de apro+imaci)n( #ue es donde se usa la f)rmula pero $omando las pare0as con$in/as a la apro+imaci)n del nivel an$erior( & #ue
corresponden cuando
h2=
h1 2
%
Trabajo de Investigación ~Métodos Numéricos~
Página 4
Los niveles de apro+imaci)n dependen de las apro+imaciones #ue se 'icieron en el nivel 3% En !eneral( si en el primer nivel( iniciamos con n apro+imaciones( en$onces alcan,aremos a lle!ar 'as$a el nivel de apro+imaci)n n%
Ejem!o" Usar el al!ori$mo de Romer! para apro+imar la in$e!ral. 2
∫e
x
ln dx
1
Con los si!uien$es se!men$os de lon!i$ud% 1 2
1 4
h1=1, h2= ,h 3= , h4 =
1 8
Trabajo de Investigación ~Métodos Numéricos~
Página 5
So!#ción" 1rimero se usa la re!la del $rapecio para llenar el nivel 4% Tenemos en$onces #ue. l ( h1 )=
2−1 1 2 [ e ln 1+ e ln 2 ] =2.560851701 2
l ( h2 )=
2−1 1 15 2 [ e ln 1+ 2 e ln1.5 + e ln 2 ] =2.189010122 4
[
(
[
(
5
3
7
2−1 1 5 3 7 l ( h3 )= e ln 1 + 2 e 4 ln + e 2 ln + e 4 ln 8 4 2 4
1
2
7
)
]
+ e2 ln 2 =2.09430857
2−1 1 1 2 7 l ( h3 )= e ln 1 + 2 e 8 ln + e 8 ln + … + e 8 ln 16 8 8 8
)
]
+ e 2 ln 2 =2.070524501
Usando las f)rmulas de Romer! para cada nivel( o$enemos la si!uien$e $ala.
2
Donde se conclu&e #ue la apro+imaci)n uscada es.
∫ e lndx ≈ 2.062586821 x
1
C#adrat#ra Ga#$$iana Es$e se convier$e en un e+celen$e m*$odo num*rico para evaluar in$e!rales definidas de funciones( por medio de suma$orias simples & f-ciles de implemen$ar% 1or o$ra par$e( es una aplicaci)n as$an$e in$eresan$e de los polinomios or$o!onales% Trabajo de Investigación ~Métodos Numéricos~
Página
Una cuadra$ura de 5auss n6pun$os llamada as" por Carl 5auss( es una cuadra$ura #ue selecciona los pun$os de la evaluaci)n de manera )p$ima & no en forma i!ualmen$e espaciada( cons$ruida para dar el resul$ado de un polinomio de !rado 7n64 o menos( ele!iles para los pun$os + & los coeficien$es
wi
para
i=1, …,n.
b
n
∫
I = f ( x ) dx =
wif ( xi ) ∑ = i 0
a
Tal cuadra$ura dar- resul$ados precisos solo si
f ( x )
es apro+imado por un
polinomio den$ro del ran!o 864(49% Si la funci)n puede ser escri$a como f ( x )=W ( x ) g ( x )
es un polinomio apro+imado & w ( x ) es conocido%
1
1
n
−1
−1
i 0
w g( x ) ∫ f ( x ) dx =∫ W ( x ) g ( x ) dx ≈ ∑ = i
i
Deri%ación de !a &órm#!a de Ga#$$'(egendre ba$ada en do$ #nto$ La cuadra$ura !aussiana de$ermina los coeficien$es de una ecuaci)n de la forma I ≅ w1 f ( x 1 ) + w 2 f ( x2 )
Donde
w1, w2
son los coeficien$es inc)!ni$as( & en con$ras$e a la re!la
$rape,oidal #ue usa pun$os e+$remos a& ( los ar!umen$os de la funci)n x 1 & x 2
a'ora no es$-n fi0os a los pun$os e+$remos( sino #ue son inc)!ni$as( & se
re#uieren de cua$ro condiciones para de$erminarlas e+ac$amen$e% Las condiciones se o$ienen suponiendo #ue la ecuaci)n de coeficien$es inde$erminados a0us$a e+ac$amen$e la in$e!ral de una cons$an$e( de una funci)n Trabajo de Investigación ~Métodos Numéricos~
Página !
lineal( e+$endiendo es$e ra,onamien$o al suponer #ue a0us$a la in$e!ral a una funci)n para)lica & a una funci)n cuica% 1
∫ 1 dx=2
w 1 f ( x 1 ) + w2 f ( x 2) =
−1
1
∫
w 1 f ( x 1 ) + w2 f ( x 2) = xdx =0 −1
1
∫
2
w 1 f ( x 1 ) + w2 f ( x 2) = x dx = −1
2 3
1
∫
3
w 1 f ( x 1 ) + w2 f ( x 2) = x dx =0 −1
Es$as ecuaciones se resuelven simul$-neamen$e( w 1 = w 2= 1
x 1=
x 2=
−1
√ 3 1
√ 3
=−0.577350269 …
=0.577350269 …
Sus$i$u&endo en la ecuaci)n de par$ida se o$iene la f)rmula de 5auss6Le!endre de los dos pun$os I ≅ w1 f ( x 1 ) + w 2 ( x 2 )= f
( √ ) ( √ ) −1 3
+ f
1
3
La f)rmula de in$e!raci)n #ue se deriv) es el miemro m-s simple de las cuadra$uras de 5auss% 1or lo $an$o( se lle!a al resul$ado in$eresan$e de #ue la
Trabajo de Investigación ~Métodos Numéricos~
Página "
suma simple de los valores de la funci)n en
( √ ) ( √ )
x =
−1 y 1 3
3
lleva una
es$imaci)n de la in$e!ral con una e+ac$i$ud de $ercer !rado%
)*todo de E#!er Es$e m*$odo se aplica para encon$rar la soluci)n a ecuaciones diferenciales ordinarias( es$o es( cuando la funci)n involucra solo una variale independien$e. dy =f ( x , y ) dx
El m*$odo se asa de forma !eneral en la pendien$e es$imada de la funci)n para e+$rapolar desde un valor an$erior a un nuevo valor. Nuevo valor =Valor anterior − Pendiente× tamañode a!o
Como $ami*n( y i+1= y1 + ∅ h
(1)
De es$a manera la formula se aplica paso a paso para encon$rar un valor en el fu$uro & as" $ra,ar la $ra&ec$oria de la soluci)n%
Trabajo de Investigación ~Métodos Numéricos~
Página #
Figura 1. Predicción de un nuevo valor en la solución. Fuente: www.gridmorelos.uaem.mx/~mcruz/cursos/mn/euler.pdf
El m*$odo de Euler u$ili,a la pendien$e al inicio del in$ervalo como una apro+imaci)n de la pendien$e promedio sore $odo el in$ervalo% La primera derivada proporciona una es$imaci)n direc$a de la pendien$e x i % = f ( x , y )
∅
f ( x , y ) ( es la ecuaci)n diferencial evaluada en x 1 &
y i
% Sus$i$u&endo es$a
es$imaci)n de la pendien$e en la ecuaci)n :4;( se $iene. y i+1= yi + f ( xi , y i ) h
(2 )
La ecuaci)n :7;( se le conoce como el m*$odo de Euler% En es$a f)rmula se predice un nuevo valor de
y
por medio de la pendien$e #ue es i!ual a la
primera derivada en el valor ori!inal de
x ( es$e nuevo valor 'ar- de
e+$rapolarse en forma lineal sore el $ama
Ejem!o" Dada la si!uien$e ecuaci)n diferencial con la condici)n inicial. y " =2 xy
y ( 0 )=1
#roximar y ( 0.5 )
Trabajo de Investigación ~Métodos Numéricos~
Página 1$
NOTA" 1rimero oservamos #ue es$a ecuaci)n s" puede resolverse por m*$odos $radicionales de ecuaciones diferenciales% 1or e0emplo( podemos aplicar el m*$odo de separaci)n de variales% =eamos las dos soluciones%
Solución Analítica+ dy =2 xy dx
dy =2 xdx y
=∫ 2 xdx ∫ dy y ln| y|= x
2
+$
Sus$i$u&endo la condici)n inicial. x =0 % y =1
ln 1=0
2
+$
0 =$
1or lo $an$o( $enemos #ue la curva soluci)n real es$- dada. lny= x
2
Trabajo de Investigación ~Métodos Numéricos~
Página 11
lny
2
x
e =e
2
x
y = e
2 por lo $an$o( el valor real #ue se pide es. y ( 0.5 )=e
(0.5 )2
=1.28403
Aplicamos el m*$odo de Euler & para ello( oservamos #ue la dis$ancia en$re x 0=0
&
x 1=0.5
no es lo suficien$emen$e pe#ue
dis$ancia en$re cinco o$enemos un valor de h =0.1 & por lo $an$o( o$endremos la apro+imaci)n deseada en cinco pasos% De es$a forma( $enemos los si!uien$es da$os.
{
x0 =0 y 0 =1 h=0,1
f ( x , y ) =2 xy
Sus$i$u&endo es$os da$os en la f)rmula de Euler( $enemos( en un primer paso.
{
x2= x0 + h =0.1
y 2= y 1+ hf ( x1 , y 1) =1 + 0.1 [ 2 ( 0.1 )( 1)] =1.02
Aplicando nuevamen$e la f)rmula de Euler( $enemos( en un se!undo paso.
{
x2 = x1 + h =0.2
y 2= y 1+ hf ( x1 , y 1) =1 + 0.1 [ 2 ( 0.1 )( 1)] =1.02
Trabajo de Investigación ~Métodos Numéricos~
Página 12
2 as" sucesivamen$e 'as$a o$ener
% Resumimos los resul$ados en la si!uien$e
$ala.
n
x n
yn
3
3
4
4
3%4
4
7
3%7
4%37
>
3%>
4%3?3@
3%
4%47B
B
3%B
4%74
Concluimos #ue el valor apro+imado( usando el m*$odo de Euler es. y ( 0.5 ) ≈ 1.2144
1ues$o #ue en es$e caso( conocemos el valor verdadero( podemos usarlo para calcular el error rela$ivo porcen$ual #ue se come$i) al aplicar la f)rmula de Euler% Tenemos #ue.
|
|∈v|=
1.28402−1.2144 × 100 1.28402
|=
5.42
,$e#docódigo" Lea a((n( Trabajo de Investigación ~Métodos Numéricos~
Página 13
' :6a;n t 0=a
y 0= &
1ara i desde 4 'as$a n 'a!a t i =a + i∗h
y i= y i−1 + h∗f ( t i−1 , y i−1 )
Fin para 1ara i desde 4 'as$a n 'a!a Imprima ( t i , y i )
-in
Siendo :$; & :&;( vec$ores( & siendo f:; un m*$odo #ue recie dos par-me$ros & arro0a un resul$ado :la funci)n;
)*todo de .e#n Es un m*$odo #ue me0ora la es$imaci)n de Euler( al es$imar la pendien$e con dos derivadas para el in$ervalo ' evaluado( una en el pun$o inicial & la o$ra en el pun$o final%
Trabajo de Investigación ~Métodos Numéricos~
Página 14
Figura 2. Corrección de la pendiente con el método de Heun al usar dos derivadas. n a! Predictor " #! corrector. Fuente: www.gridmorelos.uaem.mx/~mcruz/cursos/mn/euler.pdf
En el m*$odo de Euler la pendien$e al inicio de un in$ervalo es. y " = f ( x i+1 , y i +1)
1ara e+$rapolar linealmen$e a
y i+1
y i+ 1= yi + f ( xi , y i ) h
En el m*$odo de eun la y i+1
es una predicci)n in$ermedia conocida como
ecuaci)n predic$or% Es$a permi$e la es$imaci)n de la pendien$e al final del in$ervalo. y i+1= f ( x i+1 , y i + 1 )
Trabajo de Investigación ~Métodos Numéricos~
Página 15
Las dos pendien$es calculadas( al inicio & al final del in$ervalo se pueden cominar para o$ener una pendien$e promedio para el in$ervalo. y´ =
y + y i+ 1 2
=
f ( x i , y i ) + f ( xi +1 , y i+1 ) 2
Es$- pendien$e promedio se u$ili,a despu*s para e+$rapolar linealmen$e desde y i
( 'as$a
y i+ 1
usando el m*$odo de Euler( #ue se conoce a'ora como
ecuaci)n correc$or. y i+ 1= yi +
f ( x i , yi ) + f ( x i+1 , y i + 1 ) 2
h
La ecuaci)n posee en amos lados del si!no i!ual a y i+1 ( por lo #ue se puede aplicar de forma i$era$iva ella misma( varias veces en cada in$ervalo% De es$a manera en cada in$ervalo se podr- me0orar repe$idamen$e una es$imaci)n de y i+ 1
%
,$e#docódigo" lea f ,x 0 , y 0 , h , n
imrima x0 , y 0
arai de!de 1 ha!ta n
' 1= h∗ f ( x 0 , y 0 ) ' 2= h∗f ( x 0+ h , y 0 + ' 1 )
Trabajo de Investigación ~Métodos Numéricos~
Página 1
y 1= y 0 +
' 1 + ' 2 2
y 0= y 1 x 0= x 0+ h
imrima x0 , y 0
)*todo de Ta/!or El m*$odo de Euler no es m-s #ue un caso par$icular de los m*$odos de Ta&lor( #ue consis$en de manera !eneral en apro+imar la soluci)n por su polinomio de Ta&lor de un orden de$erminado%
{
y " =f ( x , y ) y ( x0 ) = y 0
Tal #ue presen$a /nica soluci)n
y ( x )
en un en$orno de x 0 apro+imaremos
dic'a funci)n por su polinomio de Ta&lor de orden N. ( N )
2 3 1 1 N y y ( x ) ≅ y ( x 0 ) + y " ( x 0 )( x − x0 ) + y " " ( x 0 ) ( x − x 0 ) + y " " " ( x 0 ) ( x − x 0 ) + … + x − x 0 ) ( 2 3 N (
2 el error de apro+imaci)n viene de$erminado por el res$o de orden NG4( de N + 1
manera #ue el erro proporcional a ( x − x 0 )
% Si fi0amos una sucesi)n de pun$os
e#uiespaciados. x 0 , x 1 , x 2 , … , con x n+1= x n +h ( & denominados y 0 , y 1 H a los valores apro+imados correspondien$es de y ( x ) ( $endremos #ue.
Trabajo de Investigación ~Métodos Numéricos~
Página 1!
1 1 N N 2 y n+1= y n+ y" ( x n ) h + y " " ( n ) h + … + y ( x n) h 2 N ( N + 1 Con un error de cada paso :error local 7; proporcional a h %
1ara poder aplicar el m*$odo necesi$amos conocer las derivadas de la soluci)n( pero $eniendo en cuen$a la propia ecuaci)n diferencial% y " " ( x n)= f ( x n , y n )
Ejem!o" Apli#uemos el m*$odo de Ta&lor de orden dos a la ecuaci)n y " =cos ( xy ) ( con la condici)n inicial de y ( 0 )=1 % La e+presi)n a considerar ser-. 1 2 y ( xn +1 ) ≅ y n +1= y n + y " ( x n ) h + y " " ( x n ) h 2
Tendremos en$onces #ue. x (¿ ¿ n y n ) y " " =( x n ) =f ( x n , y n ) =cos ¿
y " " =
d$o! ( xy ) =−!en ( xy ) y− !en ( xy ) y " =−!en ( xy )( y + cos ( xy )) dx
2 as" los primeros pasos $omando '3%B serian x 1= x 0+ h=0.5
2
h y 1= y 0 + h$o! ( x0 y 0 ) + (−!en ( x 0 y 0 ) ( y0+ cos ( x 0 y 0 ) ) )=¿ 2
Trabajo de Investigación ~Métodos Numéricos~
Página 1"
2
0.5 ¿ 1+ 0.5cos0 + 2
( −!en 0 ( 1 +cos ( 0 ) ) )=1.5
x 2= x 1+ h=1
2
h y 2= y 1 + h$o! ( x 1 y 1) + (−!en ( x1 y 1 ) ( y1+ cos ( x 1 y 1 ) ) )=¿ 2 2
0.5 ¿ 1.5 + 0.5cos (0.5 × 1.5 )+ 2
(−!en (0.5 × 1.5 ) ( 1.5 +cos ( 0.5 × 1.5 )) )=1.67569
2 as" sucesivamen$e% El error local en es$a apro+imaci)n ser- proporcional a
3
h
2 & por $an$o el error !loal lo ser- h %
BIB(IOGRA-0A
Trabajo de Investigación ~Métodos Numéricos~
Página 1#
Devincien,i 5% :733@;% In$e!raciones Num*ricas. In$e!raci)n de Romer!% Universidad Nacional del Nordes$e% Facul$ad de In!enier"a% 1rovincia del C'aco% Ar!en$ina% Apa,a M% :733;% 1ro!ramaci)n Aplicada% C-lculo Num*rico% Universidad Nacional de Ju0u&% Facul$ad de In!enier"a% San Salvador de Ju0u&% Ar!en$ina% Cru, M% :733;% M*$odos Num*ricos para resolver ecuaciones diferenciales Ordinarias% Universidad Au$)noma del Es$ado de Morelos% Licencia$ura en Elec$r)nica & Compu$aci)n% Cuernavaca( Morelos% M*+ico% S-nc'e, M% :733;% M*$odos Num*ricos en ecuaciones diferenciales ordinarias Tema % Universidad de Salamanca% Facul$ad de In!enier"a% Salamanca% Espa
Trabajo de Investigación ~Métodos Numéricos~
Página 2$
