{VERSION 3 0 "IBM INTEL LINUX" "3.0" } {USTYLETAB {CSTYLE "Maple Input" -1 0 "Courier" 0 1 255 0 0 1 0 1 0 0 1 0 0 0 0 }{CSTYLE "2D Math" -1 2 "Times" 0 1 0 0 0 0 0 0 2 0 0 0 0 0 0 }{CSTYLE "2D Comment" 2 18 "" 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 } {CSTYLE "" -1 256 "" 0 1 0 0 0 0 0 1 2 0 0 0 0 0 0 }{CSTYLE "" -1 257 "" 0 20 0 0 0 0 0 0 2 0 0 0 0 0 0 }{CSTYLE "" -1 258 "" 0 1 0 0 0 0 0 0 1 0 0 0 0 0 0 }{CSTYLE "" -1 259 "" 1 14 0 0 0 0 0 0 0 0 0 0 0 0 0 } {CSTYLE "" -1 260 "" 0 1 0 0 0 0 0 1 0 0 0 0 0 0 0 }{CSTYLE "" -1 261 "" 0 1 0 0 0 0 0 0 1 0 0 0 0 0 0 }{CSTYLE "" -1 262 "" 0 1 0 0 0 0 0 0 1 0 0 0 0 0 0 }{CSTYLE "" -1 263 "" 0 1 0 0 0 0 0 1 0 0 0 0 0 0 0 } {CSTYLE "" -1 264 "" 0 1 0 0 0 0 0 0 1 0 0 0 0 0 0 }{CSTYLE "" -1 265 "" 0 1 0 0 0 0 0 1 0 0 0 0 0 0 0 }{CSTYLE "" -1 266 "" 0 1 0 0 0 0 0 0 2 0 0 0 0 0 0 }{CSTYLE "" -1 267 "" 0 1 0 0 0 0 0 0 1 0 0 0 0 0 0 } {PSTYLE "Normal" -1 0 1 {CSTYLE "" -1 -1 "" 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 }0 0 0 -1 -1 -1 0 0 0 0 0 0 -1 0 }{PSTYLE "Heading 1" 0 3 1 {CSTYLE "" -1 -1 "" 1 18 0 0 0 0 0 1 0 0 0 0 0 0 0 }1 0 0 0 8 4 0 0 0 0 0 0 -1 0 }{PSTYLE "Heading 2" 3 4 1 {CSTYLE "" -1 -1 "" 1 14 0 0 0 0 0 0 0 0 0 0 0 0 0 }0 0 0 -1 8 2 0 0 0 0 0 0 -1 0 }{PSTYLE "Heading 3 " 4 5 1 {CSTYLE "" -1 -1 "" 1 12 0 0 0 0 1 0 0 0 0 0 0 0 0 }0 0 0 -1 0 0 0 0 0 0 0 0 -1 0 }{PSTYLE "Maple Plot" 0 13 1 {CSTYLE "" -1 -1 "" 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 }3 0 0 -1 -1 -1 0 0 0 0 0 0 -1 0 } {PSTYLE "" 0 256 1 {CSTYLE "" -1 -1 "" 0 1 0 0 0 0 0 0 1 0 0 0 0 0 0 } 0 0 0 -1 -1 -1 0 0 0 0 0 0 -1 0 }} {SECT 0 {PARA 3 "" 0 "" {TEXT -1 0 "" }{TEXT 256 0 "" }{TEXT 257 64 " \+ LABORATORIO DE CALCULO NUMERICO" }} {PARA 3 "" 0 "" {TEXT -1 0 "" }}{PARA 3 "" 0 "" {TEXT -1 34 " \+ " }{TEXT 258 47 "Resoluci\363n num\351rica de ecuaciones diferenciales" }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 " " 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 216 "En esta pr\341cti ca exploraremos algunos de los conceptos b\341sicos en la resoluci\363 n num\351rica de ecuaciones diferenciales, como puede ser el problema de la inestabilidad intr\355nseca de una ecuaci\363n diferencial, la \+ posible" }}{PARA 0 "" 0 "" {TEXT -1 78 "inestabilidad de un m\351todo \+ num\351rico, y el orden de convergencia de un m\351todo." }}{PARA 0 " " 0 "" {TEXT -1 0 "" }}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 49 "resta rt:with(plots): with(DEtools): with(linalg):" }}}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{SECT 1 {PARA 4 "" 0 "" {TEXT -1 49 " Rutinas auxiliar es a usar durante esta pr\341ctica." }}{PARA 0 "" 0 "" {TEXT -1 0 "" } }{PARA 0 "" 0 "" {TEXT -1 100 " Maple, por sus caracteristicas de soft ware matem\341tico eminentemente simb\363lico, tiende a enmascarar " } }{PARA 0 "" 0 "" {TEXT -1 104 "el caracter num\351rico de algunas de s us funciones. Este es el caso de la funci\363n dsolve( ) que es la qu e" }}{PARA 0 "" 0 "" {TEXT -1 105 "resuelve ecuaciones diferenciales. \+ Dada la naturaleza de nuestro curso nos interesa dejar bien claro que \+ " }}{PARA 0 "" 0 "" {TEXT -1 96 "la salida de un m\351todo num\351rico de resoluci\363n de ecuaciones diferenciales debe ser un mont\363n de " }}{PARA 0 "" 0 "" {TEXT -1 95 "numeritos, los valores y_n que el m \351todo cree que son una buena aproximaci\363n a la verdadera ( y" }} {PARA 0 "" 0 "" {TEXT -1 42 "desconocida) soluci\363n en los puntos t_ n. " }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 101 " Por ello hemos encapsulado la rutina dsolve( ) de Maple dentro de otra , diffnum( ), de forma que los " }}{PARA 0 "" 0 "" {TEXT -1 81 "par \341metros que recibe y devuelve son muy similares a la notaci\363n us ada en clase. " }}{PARA 0 "" 0 "" {TEXT -1 106 "El nuevo procedimiento recibe una ecuacion diferencial ecu, un intervalo [t0 tf], un valor \+ inicial y(t0) " }}{PARA 0 "" 0 "" {TEXT -1 95 "y un paso h, as\355 com o el nombre de un m\351todo a elegir. Los posibles m\351todos son los \+ mismos que " }}{PARA 0 "" 0 "" {TEXT -1 40 "admite dsolve( ) en su opc i\363n classical:" }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 73 " foreuler, heunform, impoly, rk2, rk3, rk4, adambash , abmoulton, etc" }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 99 "El procedimiento devuelve una matriz cuya primera fila so n los tiempos tn=t0+n*h, y la segunda los " }}{PARA 0 "" 0 "" {TEXT -1 47 "valores y_n encontrados con el m\351todo num\351rico." } {MPLTEXT 1 0 0 "" }}{PARA 0 "" 0 "" {TEXT -1 1 " " }{MPLTEXT 1 0 0 "" }{TEXT -1 0 "" }{MPLTEXT 1 0 0 "" }{TEXT -1 0 "" }}{EXCHG {PARA 0 "> \+ " 0 "" {MPLTEXT 1 0 45 "# Rutina que encapsula el interface de Maple \+ " }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 47 "# con respecto a las ecuaciones diferenciales. " }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}{PARA 0 "> \+ " 0 "" {MPLTEXT 1 0 36 "diffnum:=proc(ecu,t0,tf,y0,h,metodo)" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 1 " " }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 30 " \+ local sol,k,i,N,tt,ttl,s,yy;" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 2 " " }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 22 " N:=round((tf-t0)/h);" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 28 " ttl:=[seq(t0+k*h,k=0..N)];" }}{PARA 0 " > " 0 "" {MPLTEXT 1 0 17 " tt:=array(ttl);" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 105 " s:=dsolve(\{ecu,y(t0)=y0\},\{y(t)\},type=numeric,m ethod=classical[metodo], value=tt,stepsize=h);" }}{PARA 0 " > " 0 "" {MPLTEXT 1 0 32 " yy:=seq(s[2,1][k,2],k=1..N+1);" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 2 " " }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 37 " s ol:=array(1..2,1..N+1,[ttl,[yy]]);" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 20 " RETURN(eval(sol));" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 4 "end:" }} }{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 93 "Al devol ver los resultados de una forma no standard, podemos tener problemas p ara plotearlos." }}{PARA 0 "" 0 "" {TEXT -1 96 "La siguiente rutina re cibe una matriz generada por la rutina anterior y devuelve una estruct ura " }}{PARA 0 "" 0 "" {TEXT -1 47 "que se puede mandar directamente \+ a un plot( )." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 58 "# Rutina que recibe una matriz de dos filas y pr epara para" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 53 "# plotear la fila fil a en funci\363n de la primera (x) " }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 21 "graf := proc (A,fila)" }} {PARA 0 "> " 0 "" {MPLTEXT 1 0 18 " local data,k,N;" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 16 " N:=coldim(A);" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 49 " data := [seq([A[1,k], A[fila,k]],k = 1 .. N)];" }}{PARA 0 " > " 0 "" {MPLTEXT 1 0 16 " RETURN(data);" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 5 " end:" }}}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 97 "Finalmente en los casos en que conozcamos la soluci \363n exacta querremos compararla con la obtenida" }}{PARA 0 "" 0 "" {TEXT -1 102 "num\351ricamente. La siguiente rutina recibe una matriz \+ generada por diffnum( ) (formada por tn y yn) y " }}{PARA 0 "" 0 "" {TEXT -1 111 "una funci\363n f( ), la soluci\363n exacta. La rutina ca lcula en los diferentes tn la discrepancia entre f(tn) e yn," }} {PARA 0 "" 0 "" {TEXT -1 26 "devolviendonos un gr\341fico." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 64 "# Re cibimos una soluci\363n numerica (matriz a) y solucion ok (fun)" }} {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 27 "plot_error:=proc(a,fun,col)" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 21 " local t,y,g,k,data;" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 14 " t:=row(a,1);" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 14 " y:=row(a,2);" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 60 " data:=[seq([t[k],abs( y[k]-fun(t[k]))],k=1..coldim(sol))];" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 3 " " }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 40 " g:=plot(data,color=col,title=`Error`):" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 12 " RETURN(g);" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 4 "end :" }}}{PARA 0 "" 0 "" {TEXT -1 0 "" }}}{PARA 0 "" 0 "" {TEXT -1 0 "" } }{SECT 1 {PARA 4 "" 0 "" {TEXT -1 1 " " }{TEXT 259 46 "El caso de una \+ ecuaci\363n diferencial inestable." }}{PARA 0 "" 0 "" {TEXT -1 1 " " } }{PARA 0 "" 0 "" {TEXT -1 110 "Empezamos con un ejemplo que ilustra qu e si la soluci\363n de la ecuaci\363n diferencial original es inheren temente" }}{PARA 0 "" 0 "" {TEXT -1 99 "inestable, ning\372n m\351tod o num\351rico nos dar\341 soluciones aceptables en cualquier circustan cia. Sea " }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 81 " y'(t) = y(t) + cos(t) - sin(t). " }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 109 "Se verifica f\341cilmente que la soluci\363n exacta es y (t) = Aexp(t) + sin(t). Si imponemos como condici\363n inicial" }} {PARA 0 "" 0 "" {TEXT -1 78 " y(0)=0, la constante A debe anularse y \+ la soluci\363n \"ideal\" es y(t)=sin(t). " }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 107 "Sin embargo es f\341cil entend er que la exp(t) que anda por ah\355 agazapada nos va a traer problema s. En nuestro" }}{PARA 0 "" 0 "" {TEXT -1 109 "ejemplo hemos elegido \+ el \372nico valor de y(0) que anula la A. Por as\355 decirlo, mientra s nos movamos sobre la" }}{PARA 0 "" 0 "" {TEXT -1 108 "verdadera solu ci\363n estamos a salvo, pero en el momento en que nos apartemos un p oquito, ese error va a dar " }}{PARA 0 "" 0 "" {TEXT -1 107 "pie a que aparezca la exp(t) y evidentemente dicho t\351rmino terminar\341 ba rriendo a la \"verdadera\" soluci\363n." }}{PARA 0 "" 0 "" {TEXT -1 1 " " }{MPLTEXT 1 0 0 "" }}{PARA 0 "" 0 "" {TEXT -1 104 "Vamos a comprob ar que contra este tipo de situaciones inherentemente inestables no po demos hacer mucho. " }}{PARA 0 "" 0 "" {TEXT -1 1 " " }{MPLTEXT 1 0 0 "" }{TEXT -1 100 "Lo \372nico que podemos esperar es que usando un h p eque\361o y un m\351todo bueno las cosas funcionen m\341s o " }}{PARA 0 "" 0 "" {TEXT -1 24 "menos durante un tiempo." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 45 "#Ecuaci\363n diferencial y condiciones iniciales" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 44 "ines:= diff(y(t),t) = y(t) + cos(t) -sin(t):" }}{PARA 0 "> " 0 " " {MPLTEXT 1 0 19 "t0:=0.0: y0:=0.0: " }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 34 "# \"Verdadera\" soluci \363n y(t)=sin(t)" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 32 "ok:=proc(t) RE TURN(sin(t)); end:" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 50 "# Paso elegido y fin del intervalo de observa ci\363n:" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 17 "h:=0.5: tf:=5.0:" }} {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 25 "# Resolvemos usando euler" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 39 "so l:=diffnum(ines,t0,tf,y0,h,foreuler):" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 74 "g1:=plot([ok(t),graf(sol,2 )],t=0..tf,style=[line,point],color=[red,blue]):" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 29 "g2:=plot_error(sol,ok,green):" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 42 "# Pintamos solu ci\363n real, num\351rica y error" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 29 "display(array(1..2,[g1,g2]));" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }{TEXT -1 0 "" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}{PARA 0 "" 0 "" {TEXT -1 1 " " }{MPLTEXT 1 0 0 "" }}{PARA 0 "" 0 "" {TEXT -1 1 " " }{MPLTEXT 1 0 0 " " }}{PARA 0 "" 0 "" {TEXT -1 101 "Los puntos azules representan la sol uci\363n num\351rica (discreta). Vemos como muy pronto se alejan de la " }}{PARA 0 "" 0 "" {TEXT -1 100 "verdadera (roja). A la derecha obse rvamos como el error cometido crece exponencialmente (\277por qu\351?) ." }}{PARA 0 "" 0 "" {TEXT -1 95 " Probad con h m\341s peque\361os. S e comprueba que el error obviamente disminuye con el paso, pero " }} {PARA 0 "" 0 "" {TEXT -1 107 "eventualmente la cosa se estropea. Adem \341s el tiempo tambien se dispara, haciendo poco pr\341ctico el proce so." }}{PARA 0 "" 0 "" {TEXT -1 1 " " }{MPLTEXT 1 0 0 "" }{TEXT -1 0 " " }}{SECT 1 {PARA 4 "" 0 "" {TEXT -1 65 " La explicaci\363n: Campo de \+ direcciones de la ecuaci\363n diferencial." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 94 " Para entender porque la soluci \363n siempre parece escaparse de la verdadera y(t)=sin(t) basta " }} {PARA 0 "" 0 "" {TEXT -1 77 "volver a pintar nuestros utiles campos de direcciones. Para nuestra ecuaci\363n:" }}{PARA 0 "" 0 "" {TEXT -1 1 " " }{MPLTEXT 1 0 0 "" }{TEXT -1 0 "" }}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 16 "# Ec diferencial" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 46 "ines:= diff(y(t),t) = y(t) + cos(t) - sin(t):" }}{PARA 0 "> " 0 " " {MPLTEXT 1 0 10 "Digits:=6:" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 29 "t0 :=0.0: tf:=6.0: y0:=0.0: " }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }} {PARA 0 "> " 0 "" {MPLTEXT 1 0 22 "# Soluci\363n \"verdadera\"" }} {PARA 0 "> " 0 "" {MPLTEXT 1 0 50 "ok:=plot(sin(t),t=t0..tf,color=gree n,thickness=2):" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 20 "# Campo direcciones " }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 70 "campo:=DEplot(ines,[y],t=t0..tf,y=-1.5..1.5,color=blu e,linecolor=red):" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}{PARA 0 "> \+ " 0 "" {MPLTEXT 1 0 18 "display(campo,ok);" }{TEXT -1 0 "" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}{PARA 0 "" 0 "" {TEXT -1 0 "" } {MPLTEXT 1 0 0 "" }}{PARA 0 "" 0 "" {TEXT -1 113 "Se observa como toda s las flechitas tienden a alejarnos de la curva sin(t). La situaci\363 n es an\341loga a ir trepando" }}{PARA 0 "" 0 "" {TEXT -1 105 "por una arista muy estrecha, donde apartarnos muy poco del camino podr\355a h acernos caer por un precipicio." }}{PARA 0 "" 0 "" {TEXT -1 109 "Podem os comprobar la inestabilidad pintando varias soluciones con condicion es iniciales muy parecidas en t=0." }}{PARA 0 "" 0 "" {TEXT -1 80 "Par a ello basta incluir el parametro ini como un argumento de la funci \363n DEplot:" }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 65 "ini:=\{[0,0.001],[0,-0.001],[0,0.01],[0,-0.01],[ 0,0.1], [0,-0.1]\}:" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 74 "campo:=DEplo t(ines,[y],t=t0..tf,y=-1.5..1.5,ini,color=blue,linecolor=red):" }} {PARA 0 "> " 0 "" {MPLTEXT 1 0 18 "display(campo,ok);" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}{PARA 13 "" 1 "" {TEXT -1 0 "" }}}{PARA 0 "" 0 "" {TEXT -1 1 " " }}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}} {PARA 0 "" 0 "" {TEXT -1 100 "Claramente se observa que condiciones in iciales muy similares nos llevan a soluciones muy distintas." }}{PARA 0 "" 0 "" {TEXT -1 104 "Esto no quiere decir que el m\351todo num\351r ico haga algo intrinsecamente \"perverso\". El m\351todo no hace m\341 s " }}{PARA 0 "" 0 "" {TEXT -1 103 "que comportarse como lo har\355a u n fen\363meno f\355sico que respondiese a una ecuaci\363n diferencial \+ inestable. " }}{PARA 0 "" 0 "" {TEXT -1 103 "Es lo que se conoce como \+ el efecto mariposa (el batido de alas de una mariposa en Brasil puede \+ terminar" }}{PARA 0 "" 0 "" {TEXT -1 105 "provocando un tif\363n en Ba ngladesh). Sistemas f\355sicos que respondan a ecuaciones diferenciale s inestables " }}{PARA 0 "" 0 "" {TEXT -1 109 "son intrinsecamente imp redecibles y contradicen el ideal mecanicista del pasado siglo de que \+ el pasado y el " }}{PARA 0 "" 0 "" {TEXT -1 108 "futuro est\341n dete minados conociendo la posici\363n y velocidad de cada particula del un iverso. Por muy alta que" }}{PARA 0 "" 0 "" {TEXT -1 109 "sea la preci si\363n de los datos iniciales un sistema como el anterior siempre ter mina desbarrando (y no sabemos" }}{PARA 0 "" 0 "" {TEXT -1 8 "cuando). " }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 0 "" }} {SECT 1 {PARA 5 "" 0 "" {TEXT -1 56 " El caso contrario: una ecuaci \363n extremadamente estable:" }}{PARA 0 "" 0 "" {TEXT -1 0 "" }} {PARA 0 "" 0 "" {TEXT -1 99 " Consideremos ahora una ecuaci\363n muy \+ similar pero de un comportamiento radicalmente distinto. Sea " }} {PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 99 " \+ y'(t) = -y(t) +cos(t) + sin(t) cuaya soluci\363n es: \+ y(t) = sin(t) +Aexp(-t)" }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 108 "A diferencia de la anterior el t\351rmino que depen de de las condiciones iniciales y/o perturbaciones va ahora " }}{PARA 0 "" 0 "" {TEXT -1 105 "con exp(-t), por lo que tiende a anularse r \341pidamente. La situaci\363n ahora es an\341loga a un camino por el " }}{PARA 0 "" 0 "" {TEXT -1 100 "fondo de un valle de escarpadas lad eras. Aunque nos desviemos un poco siempre tenderemos, dejandonos" }} {PARA 0 "" 0 "" {TEXT -1 34 "caer, a volver al fondo del valle." }} {PARA 0 "" 0 "" {TEXT -1 0 "" }}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 48 "estable:= diff(y(t),t) = -y(t) + cos(t) +sin(t):" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 10 "Digits:=6:" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 29 "t0:=0.0: tf:=6.0: y0:=0.0: " }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 26 "g 2:=plot(sin(t),t=t0..tf):" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 61 "ini:= \{[0,0],[0,0.5],[0,-0.5],[0,1],[0,-1],[0,1.5], [0,-1.5]\}:" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 70 "DEplot(estable,[y],t=t0..tf,ini,y=-1.5..1.5 ,color=blue,linecolor=red);" }{TEXT -1 0 "" }}}{EXCHG {PARA 0 "> " 0 " " {MPLTEXT 1 0 0 "" }}}{PARA 0 "" 0 "" {TEXT -1 1 " " }{MPLTEXT 1 0 0 "" }{TEXT -1 1 " " }{MPLTEXT 1 0 0 "" }}}{PARA 0 "" 0 "" {TEXT -1 1 " \+ " }{MPLTEXT 1 0 0 "" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}{PARA 0 "" 0 "" {TEXT -1 1 " " }{MPLTEXT 1 0 0 "" }{TEXT -1 1 " " }{MPLTEXT 1 0 0 "" }}}{PARA 0 "" 0 "" {TEXT -1 0 "" }{MPLTEXT 1 0 0 "" }}{SECT 1 {PARA 4 "" 0 "" {TEXT -1 30 " M\351todos estables e inestables" }} {PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 108 "El probl ema de la inestabilidad ilustrado antes se refer\355a a la ecuaci\363n diferencial en s\355 misma. Veremos a " }}{PARA 0 "" 0 "" {TEXT -1 90 "continuaci\363n como un problema muy similar puede aparecer con el m\351todo num\351rico a emplear." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }} {PARA 0 "" 0 "" {TEXT -1 104 "Vamos a ilustrar a continuaci\363n la di ferencia entre un m\351todo inestable y otro estable. Para ello vamos \+ " }}{PARA 0 "" 0 "" {TEXT -1 104 "a plantear dos m\351todos num\351ric os que tienen la misma soluci\363n matem\341tica y a estudiar su compo rtamiento " }}{PARA 0 "" 0 "" {TEXT -1 49 "trabajando con precisi\363n finita. Los m\351todos son:" }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 44 "Metodo A de un s\363lo paso : y(n) = \+ " }{XPPEDIT 18 0 "y(n-1)/B;" "6#*&-%\"yG6#,&%\"nG\"\"\"\"\"\"!\"\"F)% \"BGF+" }{TEXT -1 18 " con y(0) = 1" }}{PARA 0 "" 0 "" {TEXT -1 45 "M\351todo B de dos pasos : y(n) = (" }{XPPEDIT 18 0 " B+1/B;" "6#,&%\"BG\"\"\"*&\"\"\"F%F$!\"\"F%" }{TEXT -1 50 " ) y(n-2) \+ - y(n-1) con y(0) = 1, y(1) = " }{XPPEDIT 18 0 "1/B;" "6#*&\" \"\"\"\"\"%\"BG!\"\"" }{TEXT -1 2 " ." }}{PARA 0 "" 0 "" {TEXT -1 0 " " }}{PARA 0 "" 0 "" {TEXT -1 107 " Se puede comprobar facilmente que, para esos valores iniciales, ambos m\351todos generan la misma soluc i\363n," }}{PARA 0 "" 0 "" {TEXT -1 47 "las sucesivas potencias invers as de B, esto es," }}{PARA 0 "" 0 "" {TEXT -1 75 " \+ " }{XPPEDIT 18 0 "y(n) = (1/B)^n;" "6#/-%\"yG6#%\"nG)*&\"\"\"\"\"\"%\"BG!\"\"F'" }} {PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 34 "Escribamo s dos peque\361as rutinas, " }{TEXT 263 11 "met1(B,n) " }{TEXT -1 3 "y " }{TEXT 260 11 "met2(B,n), " }{TEXT -1 75 "que calculan la potenc ia n-esima de 1/B usando la formula de ambos m\351todos:" }}{PARA 0 " " 0 "" {TEXT -1 1 " " }{MPLTEXT 1 0 0 "" }{TEXT -1 0 "" }}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 15 "met1:=proc(B,n)" }}{PARA 0 "> " 0 " " {MPLTEXT 1 0 18 " local y0,k,sig; " }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 11 " y0:=1.0; " }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 23 " for k from 2 to n do " }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 16 " sig:=y0/B; " }} {PARA 0 "> " 0 "" {MPLTEXT 1 0 13 " y0:=sig;" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 6 " od; " }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 14 " RETURN( sig);" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 5 " end:" }}}{PARA 0 "" 0 "" {TEXT -1 1 " " }{MPLTEXT 1 0 0 "" }}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 15 "met2:=proc(B,n)" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 28 " local \+ y0,y1,k,sig,factor; " }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 38 " y0:=1.0; \+ y1:=1.0/B; factor:= B+1/B;" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 23 " fo r k from 2 to n do " }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 24 " sig:=y1 *factor-y0; " }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 21 " y0:=y1; y1:=si g;" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 6 " od; " }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 14 " RETURN(sig);" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 5 " \+ end:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}{PARA 0 "" 0 " " {TEXT -1 0 "" }{MPLTEXT 1 0 1 " " }{TEXT -1 1 " " }{MPLTEXT 1 0 0 " " }{TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 102 "Comprobemos el funcio namiento de ambos m\351todos. Para ello, generaremos los 40 primeros t erminos de la " }}{PARA 0 "" 0 "" {TEXT -1 111 "soluci\363n seg\372n u no u otro m\351todo para un cierto B usando diversas precisiones (dist intos colores en la gr\341fica," }}{PARA 0 "" 0 "" {TEXT -1 118 "indic adas en la lista prec[ ] ). Los resultados se representan en una gr \341fica logaritmica en el eje de las Y, de forma" }}{PARA 0 "" 0 "" {TEXT -1 22 "que la curva y(n) = " }{XPPEDIT 18 0 "(1/B)^n;" "6#)*& \"\"\"\"\"\"%\"BG!\"\"%\"nG" }{TEXT -1 90 " debe aparecer como una r ecta log(y) = -log(B) n , de pendiente negativa (-logB.) " }} {PARA 0 "" 0 "" {TEXT -1 70 "Probemos en primer lugar con B=2 y prec ision = 4, 8, 10 digitos. " }}{PARA 0 "" 0 "" {TEXT -1 0 "" } {MPLTEXT 1 0 2 " " }}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 18 "B:=2: \+ hasta:=36: " }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 15 "prec:=[4,8,10]:" }} {PARA 0 "> " 0 "" {MPLTEXT 1 0 55 "g1:=[seq([seq([k,met1(B,k)],k=2..ha sta)],Digits=prec)]:" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 55 "g2:=[seq([s eq([k,met2(B,k)],k=2..hasta)],Digits=prec)]:" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 56 "g1:=logplot(g1,color=[red,blue,green],title=`Metodo A `):" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 56 "g2:=logplot(g2,color=[red,bl ue,green],title=`Metodo B`):" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 29 "dis play(array(1..2,[g1,g2]));" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" } {TEXT -1 0 "" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}{PARA 0 "" 0 "" {TEXT -1 0 "" }{MPLTEXT 1 0 0 "" }}{PARA 0 "" 0 "" {TEXT -1 102 "La gr\341fica de la izquierda representa los resultados del prim er m\351todo para las 3 precisiones usadas." }}{PARA 0 "" 0 "" {TEXT -1 90 "La de la derecha corresponde al segundo m\351todo. Correr ambo s m\351todos para B=2, 4, 7 y 10." }}{PARA 0 "" 0 "" {TEXT -1 97 "Anot ad en la siguiente tabla si creeis que el m\351todo est\341 haciendo l as cosas bien (OK) o mal (XX)" }}{PARA 0 "" 0 "" {TEXT -1 43 "en funci \363n de lo que veais en las gr\341ficas:" }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 4 " " }{TEXT 264 158 " B \+ 2 4 \+ 7 10 \+ " }{TEXT -1 14 " " }}{PARA 0 "" 0 "" {TEXT -1 0 "" } }{PARA 0 "" 0 "" {TEXT -1 17 " Metodo 1 " }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 12 " Metodo 2" }}{PARA 0 " " 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 99 "\277Qu\351 efecto \+ tiene aumentar la precisi\363n? \277Se arregla definitivamente el prob lema o s\363lo se retrasa?" }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 0 "" }{MPLTEXT 1 0 0 "" }{TEXT -1 104 "Lo que est\341 ocurriendo es que, apesar de ser matem\341ticamente equivalentes, el \+ primero de los m\351todos es " }}{PARA 0 "" 0 "" {TEXT -1 107 "estable , mientras que el segundo es inestable. Lo que est\341 pasando es algo muy similar a lo que vimos para " }}{PARA 0 "" 0 "" {TEXT -1 107 "la \+ ecuaci\363n diferencial inestable. Ocurre que en realidad el segundo m \351todo tiene en realidad dos posibles " }}{PARA 0 "" 0 "" {TEXT -1 113 "soluciones (antes de imponer las condiciones iniciales): una son las potencias inversas de B, (1/B)^n, mientras " }}{PARA 0 "" 0 "" {TEXT -1 214 "que la otra son sus potencias directas, B^n. Con las co ndiciones iniciales dadas la \372nica soluci\363n posible matem\341tic amente son las potencias inversas. Sin embargo, cualquier cambio en d ichas condiciones iniciales" }}{PARA 0 "" 0 "" {TEXT -1 109 "o en cual quiera de los y(n) que se van calculando introducen una peque\361a co mponente de la segunda soluci\363n, " }}{PARA 0 "" 0 "" {TEXT -1 51 "q ue termina predominando sobre los de la verdadera." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 88 "Sin embargo, hay un B pa ra el que el m\351todo 2 parece funcionar correctamente. \277Cu\341l e s?" }}{PARA 0 "" 0 "" {TEXT -1 104 "Podemos comprobar que para ese cas o las cosas siguen funcionando incluso para precisiones bajas o si nos " }}{PARA 0 "" 0 "" {TEXT -1 44 "vamos muy adelante (aumentando el par \341metro " }{TEXT 265 5 "hasta" }{TEXT -1 15 " ) en la serie." }} {PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 67 "\277Por q u\351 creeis que la inestabilidad no se manifiesta para dicho B?" }} {PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 78 "Si no sab eis responder a la pregunta anterior, haced el siguiente experimento: " }}{PARA 0 "" 0 "" {TEXT -1 25 "MAPLE tiene una funcion, " }{TEXT 256 9 "evalhf( )" }{TEXT 266 2 ", " }{TEXT -1 67 "que permite calcular una expresion usando directamente los recursos" }}{PARA 0 "" 0 "" {TEXT -1 98 "hardware del sistema (procesador matematico) en vez de la s rutinas del propio MAPLE. Repitamos el " }}{PARA 0 "" 0 "" {TEXT -1 96 "experimento anterior usando puenteando a MAPLE y usando directamen te el coprocesador matematico " }}{PARA 0 "" 0 "" {TEXT -1 108 "(ahora no usamos distintas precisiones pues ese es un par\341metro que solo \+ afecta a MAPLE). Al usar evalhf( ) " }}{PARA 0 "" 0 "" {TEXT -1 63 "es tamos usando la precisi\363n del coprocesador matem\341tico (fija):" } }{PARA 0 "" 0 "" {TEXT -1 0 "" }}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 18 "B:=2: hasta:=20: " }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }} {PARA 0 "> " 0 "" {MPLTEXT 1 0 44 "g1:=[seq([k,evalhf(met1(B,k))],k=2. .hasta)]:" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 44 "g2:=[seq([k,evalhf(met 2(B,k))],k=2..hasta)]:" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 43 "g1:=logpl ot(g1,title=`Metodo 1`,color=red):" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 44 "g2:=logplot(g2,color=blue,title=`Metodo 2`):" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 29 "display(array(1..2,[g1,g2]));" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}{PARA 13 "" 1 "" {TEXT -1 0 "" }}}{EXCHG {PARA 0 " > " 0 "" {MPLTEXT 1 0 0 "" }}}{PARA 0 "" 0 "" {TEXT -1 0 "" }{MPLTEXT 1 0 0 "" }}{PARA 0 "" 0 "" {TEXT -1 113 "Al igual que antes correr la \+ rutina para B=2, 4, 7 y 10 y anotad en la tabla vuestra opini\363n so bre si el m\351todo " }}{PARA 0 "" 0 "" {TEXT -1 53 "est\341 desbarran do (XX) o haciendo las cosas bien (OK):" }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 4 " " }{TEXT 256 3 " " }{TEXT 267 138 " B 2 4 \+ 7 10 \+ " }{TEXT -1 13 " " }}{PARA 0 "" 0 "" {TEXT -1 0 "" }} {PARA 0 "" 0 "" {TEXT -1 17 " Metodo 1 " }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 12 " Metodo 2" }}{PARA 0 " " 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 104 "\277Para qu\351 B 's parece funcionar bien ahora el segundo m\351todo? \277Que podemos d ecir tras este estudio sobre " }}{PARA 0 "" 0 "" {TEXT -1 94 "la repre sentaci\363n interna en coma flotante usada por Maple? \277y por el co procesador matem\341tico?" }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{SECT 1 {PARA 4 "" 0 "" {TEXT -1 20 " Orden de un m\351todo." }}{PARA 0 "" 0 "" {TEXT -1 1 " " }{MPLTEXT 1 0 0 "" }} {PARA 0 "" 0 "" {TEXT -1 106 " Una vez que estamos satisfechos de que \+ nuestro m\351todo es convergente (esto es, se aproxima a la soluci \363n" }}{PARA 0 "" 0 "" {TEXT -1 103 "verdadera al disminuir h) nos p lanteamos si existen m\351todos \"mejores\" que otros. Esta claro que \+ siempre" }}{PARA 0 "" 0 "" {TEXT -1 106 "puedo mejorar resultados baja ndo el paso h, pero eso tiene la clara contrapartida de aumentar el ti empo de" }}{PARA 0 "" 0 "" {TEXT -1 110 "calculo: pasar de un h dado a h/2 duplica el n\372mero de pasos y por el tiempo de ejecuci\363n. La idea es, como" }}{PARA 0 "" 0 "" {TEXT -1 67 "siempre, maximizar l os beneficios obtenidos por ese esfuerzo extra." }}{PARA 0 "" 0 "" {TEXT -1 1 " " }}{PARA 0 "" 0 "" {TEXT -1 101 "El orden de un m\351tod o estudia precisamente eso, qu\351 ganamos (cu\341nto se reduce el er ror cometido) al " }}{PARA 0 "" 0 "" {TEXT -1 110 "trabajar m\341s (d isminuir el h). Si probamos varios m\351todos y pintamos la evoluci \363n del error frente a la del" }}{PARA 0 "" 0 "" {TEXT -1 39 "paso h podremos decidir cual es mejor." }}{PARA 0 "" 0 "" {TEXT -1 90 "La t eor\355a nos dice que el error global cometido por un m\351todo num \351rico es proporcional a " }{XPPEDIT 18 0 "h^p;" "6#)%\"hG%\"pG" } {TEXT -1 14 " donde p es el" }}{PARA 0 "" 0 "" {TEXT -1 105 "orden del m\351todo. Un m\351todo de orden 1 indica que si duplicamos nuestro e sfuerzo (reducimos h a la mitad)" }}{PARA 0 "" 0 "" {TEXT -1 105 "dupl icamos los beneficios (reducimos el error a la mitad). En un m\351todo de orden 2 con la misma inversi\363n" }}{PARA 0 "" 0 "" {TEXT -1 111 "multiplicamos por cuatro los beneficios (reducimos el error a la cua rta parte). Por supuesto hay un precio a " }}{PARA 0 "" 0 "" {TEXT -1 95 "pagar, m\341s operaciones por paso, pero si los pasos pueden s er muchos menos nos va a compensar." }}{PARA 0 "" 0 "" {TEXT -1 0 "" } }{PARA 0 "" 0 "" {TEXT -1 106 "En este ejercicio vamos a determinar gr \341ficamente el orden de varios m\351todos \"inc\363gnita\". Para ell o vamos " }}{PARA 0 "" 0 "" {TEXT -1 103 "a correr dichos m\351todos c on la ecuaci\363n de siempre y'(t)=y(t), y(0)=1.0, cuya soluci\363n es y(t)=exp(t). " }}{PARA 0 "" 0 "" {TEXT -1 106 "Calcularemos el valor \+ que nos d\341 el m\351todo en t=1 para una bateria de diferentes h's \+ y lo guardaremos en " }}{PARA 0 "" 0 "" {TEXT -1 87 "una tabla. Poster iormente compararemos dichos resultados con el valor exacto y(1)=e . \+ " }}{PARA 0 "" 0 "" {TEXT -1 97 "Los h's los construimos a partir de u n h inicial, dividiendolo sucesivamente por un factor F: " } {XPPEDIT 18 0 "h;" "6#%\"hG" }{TEXT -1 4 " , " }{XPPEDIT 18 0 "h/F;" "6#*&%\"hG\"\"\"%\"FG!\"\"" }{TEXT -1 5 " , " }{XPPEDIT 18 0 "h/(F^2 );" "6#*&%\"hG\"\"\"*$%\"FG\"\"#!\"\"" }{TEXT -1 50 " , etc. \+ " }}{PARA 0 "" 0 "" {TEXT -1 96 "Para intentar estimar el orden del m\351todo representaremos los errores c ometidos frente a los h " }}{PARA 0 "" 0 "" {TEXT -1 64 "correspondie ntes. Como la relaci\363n entre el error y h tiende a :" }}{PARA 0 "" 0 "" {TEXT -1 50 " err = K " }{XPPEDIT 18 0 "h^p;" "6#)%\"hG%\"pG" }{TEXT -1 5 " " }}{PARA 0 " " 0 "" {TEXT -1 104 "lo mejor es representar dicha gr\341fica en ejes \+ logar\355tmicos. En efecto, tomando logar\355tmos en ambos ejes:" }} {PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 93 " \+ log( err ) = log(K) + p log ( h ) => y = a \+ + b x," }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 110 "Luego la gr\341fica logar\355tmica debe asemejarse a una recta, cuya pendiente b coincide con el orden p del m\351todo." }} {PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 108 "En el pa rrafo siguiente se implementa en Maple todo lo que hemos comentado. P ara medir la pendiente sobre " }}{PARA 0 "" 0 "" {TEXT -1 105 "la gr \341fica basta pinchar con el rat\363n en un par de puntos de la zona \+ recta y estimar la pendiente como " }{XPPEDIT 18 0 "Delta;" "6#%&Del taG" }{TEXT -1 4 "y / " }{XPPEDIT 18 0 "Delta;" "6#%&DeltaG" }{TEXT -1 3 "x ." }}{PARA 0 "" 0 "" {TEXT -1 4 " " }}{PARA 0 "" 0 "" {TEXT -1 79 "En primer lugar definimos la ecuaci\363n diferencial y cr eamos un vector de h[ ] :" }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 10 "Digits:=8:" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 61 "# Definimos la \+ ecu diferencial y las condiciones del problema" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 44 "#ecu:= diff(y(t),t)=1/(1+t*t) - 2*y(t)*y(t):" }} {PARA 0 "> " 0 "" {MPLTEXT 1 0 24 "ecu:= diff(y(t),t)=y(t):" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 26 "t0:=0.0: tf:=1.0: y0:=1.0:" }}{PARA 0 "> \+ " 0 "" {MPLTEXT 1 0 0 "" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 64 "# Creamo s una sucesi\363n de N valores de h, cada uno F veces menor" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 20 "F:=2: N:=5: h0:=1.0:" }}{PARA 0 "> " 0 " " {MPLTEXT 1 0 38 "h:=array(1..N,[seq(h0/F^k,k=0..N-1)]):" }}{PARA 0 " > " 0 "" {MPLTEXT 1 0 68 "printf(\"\\n H maximo %.5f ---> \+ H m\355nimo %.5f\",h[1],h[N]);" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 20 "errores=array(1..N):" }}{PARA 0 " > " 0 "" {MPLTEXT 1 0 18 "y_met=array(1..N):" }{TEXT -1 0 "" }}} {EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}{PARA 0 "" 0 "" {TEXT -1 0 "" }{MPLTEXT 1 0 0 "" }{TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 89 "Ahora definimos que m\351todo vamos a usar y pintamos el error com etido para diferentes h's:" }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 36 "# Definimos que m\351todo vamos a u sar " }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 17 "metodo:=foreuler:" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 18 "# P recisi\363n simple" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 25 "Digits:=6: pr intlevel:=0:" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}{PARA 0 "> " 0 " " {MPLTEXT 1 0 22 "for k from 1 to N do " }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 41 " sol:=diffnum(ecu,t0,tf,y0,h[k],metodo):" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 29 " # Valor \"num\351rico\" en t=1 " }} {PARA 0 "> " 0 "" {MPLTEXT 1 0 31 " y_met[k]:=sol[2,coldim(sol)]:" }} {PARA 0 "> " 0 "" {MPLTEXT 1 0 3 "od:" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 42 "# Valor real calculado con mucha precisi\363n" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 28 " yok:=eval f(exp(1.0),24): " }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 3 " " }}{PARA 0 " > " 0 "" {MPLTEXT 1 0 59 "# Discrepancia entre los valores encontrados y el verdadero" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 20 "for k from 1 to \+ N do" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 32 " errores[k]:=abs(yok-y_met [k]);" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 5 "od; " }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 40 "data:=[seq([h[k ],errores[k]],k = 1..N)]:" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 70 "loglog plot(data,color=red,labels=[`Intervalos: h_n`,`Errores: E_n` ]):" } {TEXT -1 0 "" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}{PARA 0 "" 0 "" {TEXT -1 1 " " }{MPLTEXT 1 0 0 "" }}{PARA 0 "" 0 "" {TEXT -1 111 "Vemos que efectivamente el resultado (sobre todo a la izquierd a, para h peque\361as) se ajusta bastante bien a una" }}{PARA 0 "" 0 " " {TEXT -1 113 "recta. Pinchando con el rat\363n hallar un par de punt os y a partir de ellos estimar la pendiente (orden) del m\351todo" }} {PARA 0 "" 0 "" {TEXT -1 115 "(Nota: para facilitar las cuentas proba d a seleccionar dos puntos separados por 1/2 en el eje de las X y ento nces " }}{PARA 0 "" 0 "" {TEXT -1 81 "la pendiente a estimar ser\341 s implemente 2 veces la variaci\363n en el eje de las Y)." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 92 " Repetir para los m \351todos incognita heunform y rk3. \277De qu\351 orden resulta ser \+ cada m\351todo?" }}{PARA 0 "" 0 "" {TEXT -1 98 "Anotad tambi\351n el o rden del error m\355nimo que se alcanza con cada m\351todo (pinchando \+ con el rat\363n y " }}{PARA 0 "" 0 "" {TEXT -1 76 "anotando el valor e n las y's, recordando que el error sera 10 a dicho valor)" }}{PARA 0 " " 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 92 " \+ x1 y1 x2 y2 \+ " }{XPPEDIT 18 0 "Delta;" "6#%&DeltaG" }{TEXT -1 11 "x " }{XPPEDIT 18 0 "Delta;" "6#%&DeltaG" }{TEXT -1 15 "y m= " } {XPPEDIT 18 0 "Delta;" "6#%&DeltaG" }{TEXT -1 4 "y / " }{XPPEDIT 18 0 "Delta;" "6#%&DeltaG" }{TEXT -1 34 "x p M\355n Er ror" }}{PARA 0 "" 0 "" {TEXT -1 87 "__________________________________ _____________________________________________________" }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 143 "foreuler -> \+ -1.20 -1.09 -0.70 -0.64 0.5 0.45 \+ 0.9 1 ? 1e-1.5" }}{PARA 0 "" 0 " " {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 13 "heunform -> " }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 16 "rk3 - >" }}{PARA 0 "" 0 "" {TEXT -1 0 "" }{MPLTEXT 1 0 0 "" }}{PARA 0 "" 0 " " {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 5 "" 0 "" {TEXT -1 40 "El problema de hacer h demasiado peque\361a" }}{PARA 0 " " 0 "" {TEXT -1 2 " " }}{PARA 0 "" 0 "" {TEXT -1 103 "Como se ha vist o en los ejemplos anteriores el error obviamente depende del h escogid o. Una vez que nos" }}{PARA 0 "" 0 "" {TEXT -1 101 "decidamos por un m \351todo (siempre que fuese convergente) podr\355amos pensar que como \+ el error tiende a " }}{PARA 0 "" 0 "" {TEXT -1 9 "cero con " } {XPPEDIT 18 0 "h;" "6#%\"hG" }{TEXT -1 89 " podremos conseguir una pr ecisi\363n arbitrariamente grande haciendo el paso m\341s peque\361o. \+ " }}{PARA 0 "" 0 "" {TEXT -1 110 "Veamos si es as\355. Para ello bas ta aumentar la constante N que determinaba el numero de h's a usar; \+ podemos" }}{PARA 0 "" 0 "" {TEXT -1 62 "pasar de N=5 a N=10, usando \+ un h m\355nimo 32 veces mas peque\361o." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 91 "Correr el programa con ese nuevo N ( esto es un h_min de aproximadamente 0.002) y anotad el " }}{PARA 0 "" 0 "" {TEXT -1 52 "error m\355nimo alcanzado para los siguientes m\351t odos: " }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT 262 106 "METODO: foreuler heunform \+ rk3 rk4" }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 11 "ERR MIN: " }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 110 "\277Se observa alg\372n desviamiento de lo esperado? \277Existe alg\372n tipo de barrera que \+ ning\372n m\351todo puede superar?" }}{PARA 0 "" 0 "" {TEXT -1 57 " \277Por qu\351 los m\351todos de orden alto parecen desvariar m\341s? " }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 107 "Como parece que rk4 es quien nos proporciona el error m\341s bajo, mantene dlo como m\351todo y variar el n\372mero " }}{PARA 0 "" 0 "" {TEXT -1 60 "de d\355gitos (precisi\363n de la m\341quina). Anotad los resul tados:" }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 25 "METODO= rk4 " }{TEXT 261 92 "DIGITS \+ 4 8 12 16 " }} {PARA 0 "" 0 "" {TEXT -1 8 " " }}{PARA 0 "" 0 "" {TEXT -1 52 " \+ ERR MIN " }}{PARA 0 "" 0 " " {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 84 "\277Qu\351 podemos conc luir de estos resultados? \277Cu\341l es el l\355mite pr\341ctico del \+ h a usar?" }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{SECT 1 {PARA 4 "" 0 "" {TEXT -1 39 " Sumario: comparaci\363n de varios m\351todos" }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 113 "En este apartado simplemente se trata de resumir lo estudiado en el punto anterior jun tando en una \372nica gr\341fica " }}{PARA 0 "" 0 "" {TEXT -1 90 "el c omportamiento de varios m\351todos. Partimos de la misma ecuaci\363n d iferencial de siempre:" }}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 50 "ec u:= diff(y(t),t)=y(t): #1/(1+t*t) -2*y(t)*y(t):" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 26 "t0:=0.0: tf:=1.0: y0:=1.0:" }}}{PARA 0 "" 0 "" {TEXT -1 0 "" }{MPLTEXT 1 0 0 "" }{TEXT -1 0 "" }{MPLTEXT 1 0 0 "" }}{PARA 0 "" 0 "" {TEXT -1 108 "A continuaci\363n generamos un gr\341fico como el anterior, pero para varios m\351todos a la vez, y con una precisi \363n" }}{PARA 0 "" 0 "" {TEXT -1 14 "de 18 digitos:" }}{PARA 0 "" 0 " " {TEXT -1 0 "" }}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 37 "N_metodos: =5: g:=array(1..N_metodos):" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 36 "# En umeramos los m\351todos a comparar " }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 37 "met:=\{foreuler,rk2,rk3,rk4, adambash\}:" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 34 "col:=\{red,green,blu e,cyan,yellow\}:" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}{PARA 0 "> \+ " 0 "" {MPLTEXT 1 0 18 "# Precisi\363n usada:" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 11 "Digits:=18:" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }} {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 14 "printlevel:=0:" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 19 "# Creamos las h[N] " }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 7 "N:=10: " }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 54 "F:=2: h0:=1.0: h:=array(1..N,[seq(h0/F^k,k=0..N-1 )]):" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 18 "error=array(1..N):" }} {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 42 "# Valor real calculado con m ucha precisi\363n" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 48 " yok:=evalf( exp(1.0),24): #yok:=evalf(0.5,32):" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 2 " " }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 33 "# barremos los diferentes \+ m\351todos" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 28 "for m from 1 to N_met odos do" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 43 " printf(\"**** Metodo %s *****\\n\",met[m]); " }} {PARA 0 "> " 0 "" {MPLTEXT 1 0 53 " printf(\" k Paso(h) e[h] Pri mer Digito mal\\n\");" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 29 "# barremos las diferentes h's" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 24 " for k from 1 to N do " }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 41 " a:=diffnum(ecu, t0,tf,y0,h[k],met[m]):" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 1 " " }} {PARA 0 "> " 0 "" {MPLTEXT 1 0 48 " error[k]:=evalf(abs(yok-a[2,col dim(a)]),24);" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 83 " printf(\" %2d \+ %06.4f %5.2e %0.f\\n\",k,h[k],error[k],-log10(error[k])); " }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 2 " " }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 56 " data:=[seq([log10(h[k]),log10(error[k])],k = 1..N)]:" }} {PARA 0 "> " 0 "" {MPLTEXT 1 0 4 " " }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 39 " g[m]:=plot(data,color=col[m]): " }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 7 " od; " }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 3 "od;" }} {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 4 "met;" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 4 "col;" }}{PARA 0 "> " 0 " " {MPLTEXT 1 0 26 "# Pintamos los resultados:" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 34 "display(seq(g[k],k=1..N_metodos));" }{TEXT -1 0 "" }} }{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}{PARA 0 "" 0 "" {TEXT -1 0 "" }{MPLTEXT 1 0 0 "" }{TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 111 "Probad a cambiar la precisi\363n de los c\341lculos cambiando Di gits (p.e.. 14, 10, 6) y estudiar cual es el l\355mite" }}{PARA 0 "" 0 "" {TEXT -1 112 " inferior \"fiable\" del error y el comportamiento de las distintas curvas. Anotad en la tabla inferior la mejor" }} {PARA 0 "" 0 "" {TEXT -1 105 " precisi\363n (digito donde aparece el \+ error) encontrada para cada m\351todo, indicando si pensais que podr \355a " }}{PARA 0 "" 0 "" {TEXT -1 68 "seguir aumentando (bajando h) c on un OK o si se ha estancado (XX)." }}{PARA 0 "" 0 "" {TEXT -1 75 " \277El error m\355nimo alcanzado depende del m\351todo o m\341s bien d e la precisi\363n? " }}{PARA 0 "" 0 "" {TEXT -1 52 "\277Cu\341ndo sur gen anomalias al ideal de la l\355nea recta?" }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 256 "" 0 "" {TEXT -1 6 "METODO" }{TEXT 256 111 " \+ foreuler rk2 rk3 \+ rk4 adambash " }}{PARA 0 "" 0 "" {TEXT -1 85 "__ ______________________________________________________________________ _____________" }}{PARA 0 "" 0 "" {TEXT -1 134 "DIGITS 18 \+ 4 (OK) 6 (OK) 10 (OK) \+ 13 (OK) 15 (OK)" }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 29 " 14" }} {PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 29 " \+ 10" }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 " " {TEXT -1 31 " 6 " }{MPLTEXT 1 0 1 " " } {TEXT -1 0 "" }}}{PARA 0 "" 0 "" {TEXT -1 0 "" }{MPLTEXT 1 0 0 "" } {TEXT -1 0 "" }}}{PARA 0 "" 0 "" {TEXT -1 0 "" }{MPLTEXT 1 0 0 "" } {TEXT -1 0 "" }{MPLTEXT 1 0 0 "" }}}{MARK "5 0" 146 }{VIEWOPTS 1 1 0 1 1 1803 }