{VERSION 3 0 "IBM INTEL NT" "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 "2D Input" 2 19 "" 0 1 255 0 0 1 0 0 0 0 0 0 0 0 0 }{CSTYLE "2 D Output" 2 20 "" 0 1 0 0 255 1 0 0 0 0 0 0 0 0 0 }{CSTYLE "" -1 256 " " 1 14 0 0 0 0 0 1 0 0 0 0 0 0 0 }{CSTYLE "" -1 257 "" 0 1 0 0 0 0 0 2 0 0 0 0 0 0 0 }{CSTYLE "" -1 258 "" 0 1 0 0 0 0 0 0 2 0 0 0 0 0 0 } {CSTYLE "" -1 259 "" 0 1 0 0 0 0 0 0 2 0 0 0 0 0 0 }{CSTYLE "" -1 260 "" 1 12 0 0 0 0 0 2 0 0 0 0 0 0 0 }{CSTYLE "" -1 261 "" 1 10 0 0 0 0 0 1 2 0 0 0 0 0 0 }{CSTYLE "" -1 262 "" 0 10 0 0 0 0 0 2 2 0 0 0 0 0 0 }{CSTYLE "" -1 263 "" 1 12 0 0 0 0 0 1 2 0 0 0 0 0 0 }{CSTYLE "" -1 264 "" 1 14 0 0 0 0 0 1 0 0 0 0 0 0 0 }{CSTYLE "" -1 265 "" 0 1 0 0 0 0 0 1 1 0 0 0 0 0 0 }{CSTYLE "" -1 266 "" 1 18 0 0 0 0 0 0 0 0 0 0 0 0 0 }{CSTYLE "" -1 267 "" 0 20 0 0 0 0 0 0 2 0 0 0 0 0 0 }{CSTYLE "" -1 268 "" 1 14 0 0 0 0 0 0 0 0 0 0 0 0 0 }{CSTYLE "" -1 269 "" 0 1 0 0 0 0 0 1 1 0 0 0 0 0 0 }{CSTYLE "" -1 270 "" 0 1 0 0 0 0 0 2 2 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 "Text Output" -1 2 1 {CSTYLE "" -1 -1 "Courier" 1 10 0 0 255 1 0 0 0 0 0 1 3 0 3 }1 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 "" 2 6 1 {CSTYLE "" -1 -1 "" 0 1 0 0 0 0 0 0 0 0 0 0 2 0 0 }0 0 0 -1 -1 -1 0 0 0 0 0 0 -1 0 }{PSTYLE "Warning" 2 7 1 {CSTYLE "" -1 -1 "" 0 1 0 0 255 1 0 0 0 0 0 0 1 0 0 }0 0 0 -1 -1 -1 0 0 0 0 0 0 -1 0 }{PSTYLE "Maple Output" 0 11 1 {CSTYLE "" -1 -1 "" 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 }3 3 0 -1 -1 -1 0 0 0 0 0 0 -1 0 }{PSTYLE "M aple 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 1 0 0 0 0 0 0 0 }3 0 0 -1 -1 -1 0 0 0 0 0 0 -1 0 } {PSTYLE "" 4 257 1 {CSTYLE "" -1 -1 "" 1 12 0 0 0 0 0 2 1 0 0 0 0 0 0 }0 0 0 -1 -1 -1 0 0 0 0 0 0 -1 0 }{PSTYLE "" 0 258 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 9 " " }{TEXT 256 0 "" }{TEXT -1 31 "LABORATORIO DE CALCULO NUMERICO" }}{PARA 0 "" 0 "" {TEXT -1 0 " " }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 3 "" 0 "" {TEXT -1 6 " \+ " }{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 248 "En esta pr\341ctica exploraremos algunos de los conceptos b\341sicos en la resoluci\363n num\351rica de ecuaci ones diferenciales, como pueden ser la aproximaci\363n de la soluci \363n num\351rica a la soluci\363n exacta y la velocidad de esa aproxi maci\363n: orden del m\351todo." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }} {PARA 0 "" 0 "" {TEXT -1 0 "" }}{SECT 1 {PARA 4 "" 0 "" {TEXT -1 53 " \+ 0.- Rutinas auxiliares a usar durante esta pr\341ctica." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 1 " " }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{SECT 1 {PARA 3 "" 0 "" {TEXT -1 1 " " }{TEXT 264 63 "1.- Aproximaci\363n de la soluci\363n num\351rica a la soluci\363n exacta. " }}{PARA 3 "" 0 "" {TEXT 263 0 "" }{TEXT 261 0 "" }{TEXT -1 0 "" }}}{PARA 4 "" 0 "" {TEXT -1 0 "" }} {SECT 1 {PARA 4 "" 0 "" {TEXT -1 22 " 2.- Orden del m\351todo." }} {PARA 0 "" 0 "" {TEXT -1 0 "" }}}{PARA 0 "" 0 "" {TEXT -1 1 " " }} {SECT 1 {PARA 4 "" 0 "" {TEXT -1 5 " 3.- " }{TEXT 256 10 "Ejercicio." }}{PARA 4 "" 0 "" {TEXT -1 0 "" }}}{PARA 0 "" 0 "" {TEXT -1 0 "" }} {PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 309 "Otros as pectos relacionados con la resoluci\363n num\351rica de ecuaciones dif erenciales ordinarias son el problema de la estabilidad/inestabilidad \+ de los m\351todos num\351ricos, diversas t\351cnicas adaptativas y est rategias de tipo Predictor-Corrector. A continuaci\363n se presentan e jemplos ilustrativos de estos aspectos." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{SECT 1 {PARA 4 "" 0 "" {TEXT -1 35 " 4.- M\351todos estables e inestables." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 104 "Vamos a ilustrar a continuaci \363n la diferencia entre un m\351todo inestable y otro estable. Para \+ ello vamos " }}{PARA 0 "" 0 "" {TEXT -1 104 "a plantear dos m\351todos num\351ricos que tienen la misma soluci\363n matem\341tica y a estudi ar su comportamiento " }}{PARA 0 "" 0 "" {TEXT -1 49 "trabajando con p recisi\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 soluci\363n:" }}{PARA 0 "" 0 "" {TEXT -1 47 "las sucesivas p otencias inversas 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 "Escribamos dos peque\361as rutinas, " }{TEXT 258 11 "met1(B,n) \+ " }{TEXT -1 3 "y " }{TEXT 257 11 "met2(B,n), " }{TEXT -1 75 "que calc ulan la potencia n-esima de 1/B usando la formula de ambos m\351todos: " }}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}{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 1 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:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}{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 " for 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:=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:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}{PARA 0 "" 0 "" {TEXT -1 0 "" }{MPLTEXT 1 0 0 "" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 21 "met1(2,2);met2(10, 2);" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}{EXCHG {PARA 0 "> " 0 " " {MPLTEXT 1 0 0 "" }}{PARA 0 "" 0 "" {TEXT -1 1 " " }{MPLTEXT 1 0 0 " " }{TEXT -1 0 "" }}}{PARA 0 "" 0 "" {TEXT -1 102 "Comprobemos el funci onamiento de ambos m\351todos. Para ello, generaremos los 40 primeros \+ terminos de la " }}{PARA 0 "" 0 "" {TEXT -1 111 "soluci\363n seg\372n \+ uno u otro m\351todo para un cierto B usando diversas precisiones (dis tintos colores en la gr\341fica," }}{PARA 0 "" 0 "" {TEXT -1 118 "indi cadas 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. " }}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}{PARA 0 "" 0 "" {TEXT -1 0 "" }{MPLTEXT 1 0 2 " " }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 18 "B:=2: hasta:=40: " }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 15 "prec:=[6,8,10]:" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}{PARA 0 " > " 0 "" {MPLTEXT 1 0 84 "printf(\" METODO A \+ METODO B\\n\");" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 90 "printf(\" DIGITOS: 6 8 10 \+ 6 8 10\\n\");" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 29 "for n from 4 to hasta by 4 do" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 26 " printf(\"n=%3d - > \",n);" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 21 " for k from 1 to 3 d o" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 21 " Digits:=prec[k]; " }} {PARA 0 "> " 0 "" {MPLTEXT 1 0 31 " printf(\"%6.2e \",met1(B,n));" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 4 " od;" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 22 " printf(\"******** \");" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 21 " for k from 1 to 3 do" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 19 " Digits:=prec[k];" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 31 " p rintf(\"%6.2e \",met2(B,n));" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 4 " od ;" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 14 " printf(\"\\n\");" }}{PARA 0 " > " 0 "" {MPLTEXT 1 0 3 "od;" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }} {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 61 "g1:=[seq([seq([4*k,met1(B,4*k)],k=1..hasta/4)],Digits=prec)]:" }} {PARA 0 "> " 0 "" {MPLTEXT 1 0 61 "g2:=[seq([seq([4*k,met2(B,4*k)],k=1 ..hasta/4)],Digits=prec)]:" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 56 "g1:=l ogplot(g1,color=[red,blue,green],title=`Metodo A`):" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 56 "g2:=logplot(g2,color=[red,blue,green],title=`Metod o B`):" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 29 "display(array(1..2,[g1,g2 ]));" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}{PARA 0 "" 0 " " {TEXT -1 109 "Los resultados se presentan num\351rica y gr\341ficame nte. A la izquierda se presentan los resultados para las tres" }} {PARA 0 "" 0 "" {TEXT -1 92 "precisiones usando el primer m\351todo. A la derecha tenemos los resultados del segundo m\351todo." }}{PARA 0 " " 0 "" {TEXT -1 104 " Corred ambos m\351todos para B=2, 4, 7 y 10. \+ Anotad en la siguiente tabla si creeis que el m\351todo est\341 " }} {PARA 0 "" 0 "" {TEXT -1 101 "haciendo las cosas bien (OK) o mal (XX) \+ en funci\363n de lo que veais en las gr\341ficas y en los listados." } }{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 4 " " } {TEXT 259 76 " 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 " Metod o 2" }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 99 " \277Qu\351 efecto tiene aumentar la precisi\363n? \277Se arregla defin itivamente el problema 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\341ticament e equivalentes, el primero de los m\351todos es " }}{PARA 0 "" 0 "" {TEXT -1 107 "estable, mientras que el segundo es inestable. Lo que es t\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 condici ones iniciales): una son las potencias inversas de B, (1/B)^n, mientr as " }}{PARA 0 "" 0 "" {TEXT -1 108 "que la otra son sus potencias dir ectas, B^n. Con las condiciones iniciales dadas la \372nica soluci \363n posible " }}{PARA 0 "" 0 "" {TEXT -1 106 "matem\341ticamente son las potencias inversas. Sin embargo, cualquier cambio en dichas condi ciones iniciales " }}{PARA 0 "" 0 "" {TEXT -1 108 "o en cualquiera de \+ los y(n) que se van calculando introducen una peque\361a componente de la segunda soluci\363n, " }}{PARA 0 "" 0 "" {TEXT -1 52 "que termina \+ predominando sobre los de la verdadera." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 88 "Sin embargo, hay un B para el que el m\351todo 2 parece funcionar correctamente. \277Cu\341l es?" }} {PARA 0 "" 0 "" {TEXT -1 104 "Podemos comprobar que para ese caso las \+ cosas siguen funcionando incluso para precisiones bajas o si nos" }} {PARA 0 "" 0 "" {TEXT -1 44 "vamos muy adelante (aumentando el par\341 metro " }{TEXT 260 5 "hasta" }{TEXT -1 15 " ) en la serie." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 67 "\277Por qu\351 cr eeis que la inestabilidad no se manifiesta para dicho B?" }}{PARA 0 " " 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 78 "Si no sabeis respo nder a la pregunta anterior, haced el siguiente experimento:" }}{PARA 0 "" 0 "" {TEXT -1 25 "MAPLE tiene una funcion, " }{TEXT 256 9 "evalhf ( )" }{TEXT 261 2 ", " }{TEXT -1 67 "que permite calcular una expresio n usando directamente los recursos" }}{PARA 0 "" 0 "" {TEXT -1 98 "har dware del sistema (procesador matematico) en vez de las rutinas del pr opio MAPLE. Repitamos el " }}{PARA 0 "" 0 "" {TEXT -1 96 "experimento \+ anterior usando puenteando a MAPLE y usando directamente el coprocesad or matematico " }}{PARA 0 "" 0 "" {TEXT -1 108 "(ahora no usamos disti ntas precisiones pues ese es un par\341metro que solo afecta a MAPLE). Al usar evalhf( ) " }}{PARA 0 "" 0 "" {TEXT -1 63 "estamos 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: hast a:=40: " }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 13 "printf(\"\\n\");" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 64 "printf(\" METODO A METODO B\\n\") ;" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 29 "for n from 4 to hasta by 4 do" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 24 " printf(\" \");" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 27 " printf(\"n=%3d -> \",n); " }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 73 " \+ printf(\"%6.2e ****** %6.2e\\n\",evalhf(met1(B,n)),evalhf(met2(B,n) )); " }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 3 "od;" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}{PARA 0 "> \+ " 0 "" {MPLTEXT 1 0 50 "g1:=[seq([4*k,evalhf(met1(B,4*k))],k=1..hasta/ 4)]:" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 50 "g2:=[seq([4*k,evalhf(met2(B ,4*k))],k=1..hasta/4)]:" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 43 "g1:=logp lot(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 0 "" 0 "" {TEXT -1 0 "" }}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}{PARA 0 "" 0 "" {TEXT -1 113 "Al igual qu e antes correr la rutina para B=2, 4, 7 y 10 y anotad en la tabla vue stra opini\363n sobre si el m\351todo " }}{PARA 0 "" 0 "" {TEXT -1 53 "est\341 desbarrando (XX) o haciendo las cosas bien (OK):" }}{PARA 0 " " 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 4 " " }{TEXT 256 3 " " }{TEXT 262 72 " 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 " Metod o 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 decir tras este estudio sobre " }}{PARA 0 "" 0 "" {TEXT -1 94 "la representaci\363n interna en coma flotante usada por M aple? \277y por el coprocesador matem\341tico?" }}}{SECT 1 {PARA 4 "" 0 "" {TEXT -1 67 " 5.- Extrapolaci\363n de Richardson: aplicaci\363n a m\351todos adaptativos." }}{EXCHG {PARA 0 "" 0 "" {TEXT -1 0 "" }} {PARA 256 "" 0 "" {TEXT 258 0 "" }{TEXT 270 0 "" }{TEXT -1 0 "" } {TEXT 257 215 "A continuaci\363n vamos a hacer uso de la aplicaci\363n m\341s importante de la extrapolaci\363n de Richardson: utilizar la e stimaci\363n del error local que estamos cometiendo para construir un \+ m\351todo adaptativo de paso variable." }}{PARA 257 "" 0 "" {TEXT 259 400 " Dado un m\351todo cualquiera a usar, la idea es hacer cada paso \+ dos veces, una de paso h y otra de dos pasos de tama\361o h/2. Una com paraci\363n de los resultados nos dar\341 una estimaci\363n del error \+ cometido en dicho paso. Dicha estimaci\363n se compara con un umbral e stablecido por el usuario. Si es similar el paso h se mantiene. Si es \+ mucho mayor se reduce el paso, mientras que si es mucho menor se aumen ta." }}{PARA 0 "" 0 "" {TEXT -1 117 "La idea es tomar pasos diferentes , largos en zonas \"faciles\", cortos en zonas \"dif\355ciles\", repa rtiendo asi el error." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 95 "Para comparar m\351todos adaptativos vs. paso fijo vam os a usar la siguiente ecuaci\363n diferencial:" }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 32 " y' (t) = " }{XPPEDIT 18 0 "1/(1+t^2);" "6#*&\"\"\"\"\"\",&\"\"\"F%*$%\"t G\"\"#F%!\"\"" }{TEXT -1 7 " - 2 " }{XPPEDIT 18 0 "y(t)^2;" "6#*$-% \"yG6#%\"tG\"\"#" }{TEXT -1 48 " , con y(0)=0 , cuya soluci\363n ex acta es " }{XPPEDIT 18 0 "y(t);" "6#-%\"yG6#%\"tG" }{TEXT -1 4 " \+ = " }{XPPEDIT 18 0 "t/(1+t^2);" "6#*&%\"tG\"\"\",&\"\"\"F%*$F$\"\"#F% !\"\"" }{TEXT -1 5 " " }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 " " 0 "" {TEXT -1 110 "especialmente interesante porque presenta cambios muy r\341pidos de pendiente, que nos serviran para ilustrar las" }} {PARA 0 "" 0 "" {TEXT -1 92 "ventajas de un m\351todo adaptativo. Veam os en primer lugar la soluci\363n usando un Euler normal:" }}{PARA 0 " " 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 23 "# Ecuacion diferencial " }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 40 "ecu:=diff(y(t),t)=1/(1+t*t)-2*y(t)*y(t):" }} {PARA 0 "> " 0 "" {MPLTEXT 1 0 20 "# Soluci\363n verdadera" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 35 "ok:=proc(t) RETURN(t/(1+t*t)); end:" }} {PARA 0 "> " 0 "" {MPLTEXT 1 0 35 "# Condiciones iniciales e intervalo " }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 21 "t0:=0:tf:=10:y0:=0.0:" }}} {PARA 0 "" 0 "" {TEXT -1 0 "" }}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 7 "h:=0.5:" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 38 "eul:=diffnum(ecu,t0,t f,y0,h,foreuler):" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 72 "plot([ok(t),gr af(eul,2)],t=t0..tf,style=[line,point],color=[red,blue]);\n" }}}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" } }}{PARA 0 "" 0 "" {TEXT -1 103 "Se comprueba que el intervalo escogido h=0.5 es demasiado grande y no somos capaces de seguir el brusco" }} {PARA 0 "" 0 "" {TEXT -1 101 "cambio de la funci\363n en t=1. Por el c ontrario, el ajuste para t's mayores de 2 es bastante razonable." }} {PARA 0 "" 0 "" {TEXT -1 109 "Probad con valores m\341s peque\361os de h (0.1, 0.05, etc.). \277Para qu\351 valores de h empieza a ser bueno el ajuste " }}{PARA 0 "" 0 "" {TEXT -1 81 "en la c\372spide, alrededo r de t=1? \277Qu\351 pasa entonces para la zona de la derecha? " }} {PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 85 "El ejempl o anterior ilustra una situaci\363n muy com\372n. T\355picamente una s oluci\363n exhibe " }{TEXT 256 12 "transitorios" }{TEXT -1 2 ", " }} {PARA 0 "" 0 "" {TEXT -1 101 "comportamientos muy localizados que desa parecen con el tiempo. Para que un m\351todo de paso fijo \"vea\" " }} {PARA 0 "" 0 "" {TEXT -1 102 "bien dichos transitorios debemos usar un paso muy peque\361o, paso que ser\341 innecesariamente peque\361o una " }}{PARA 0 "" 0 "" {TEXT -1 103 "vez pasado el transitorio. As\355 s urge la idea de implementar m\351todos que adapten su paso a la soluci \363n. " }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 108 "En la vida real, un m\351todo adaptativo tiene muchas sutilezas, \+ acerca de la regla de actualizaci\363n del paso, " }}{PARA 0 "" 0 "" {TEXT -1 102 "posibles valores m\341ximos y m\355nimos para el paso, e tc. Nosotros vamos a programar un m\351todo muy simple," }}{PARA 0 "" 0 "" {TEXT -1 211 "basado en un Euler. Aunque poco eficaz en la pr\341 ctica por su bajo orden y poca sofisticaci\363n, va a servir para ilus trar las ideas b\341sicas detr\341s de un m\351todo adaptativo, compar andolo con el original sin adaptar." }}{PARA 0 "" 0 "" {TEXT -1 0 "" } }{PARA 0 "" 0 "" {TEXT -1 373 "La siguiente rutina implementa un adapt ativo basado en aplicar Euler , primero con unico paso h y luego con d os pasos de tama\361o h/2. A partir de la discrepancia (extrapolaci \363n de Richardson) decido autom\341ticamente si aumentar o disminuir el paso. Un control adicional me marca que nunca disminuya el paso po r un factor mayor que 8 ni lo aumente por un factor mayor que 2. " }} {PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 111 "Adem\341 s de los par\341metros usuales le suministraremos un paso inicial h0 \+ y un umbral frente al cual comparar la " }}{PARA 0 "" 0 "" {TEXT -1 21 "estimaci\363n del error:" }}{PARA 0 "" 0 "" {TEXT -1 1 " " } {MPLTEXT 1 0 0 "" }}{PARA 0 "" 0 "" {TEXT -1 103 "En primer lugar defi nimos la f(x,y) de nuestra ecuaci\363n diferencial, usada por nuestra \+ rutina. Despues " }}{PARA 0 "" 0 "" {TEXT -1 43 "escribimos el c\363di go del m\351todo adaptativo:" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" } {TEXT -1 0 "" }{MPLTEXT 1 0 0 "" }{TEXT -1 0 "" }}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 8 "# f(x,y)" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 50 "de riv := proc (x,y) RETURN(1/(1+x*x)-2*y*y); end:" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 25 "# Rutina Euler \+ adaptativo" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 36 "euler_adapta:= proc(t 0,tf,yini,h0,U)" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 74 "local t,k,h,ylas t,ynext,y1,m1,m2,err,sol,Np,h2,y2,ytemp,F,hnext,Uh,abserr;" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 25 "sol:= array(1..4,1..1000);" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 56 "sol[1,1]:=t 0; sol[2,1]:=y0; sol[3,1]:=h0; sol[4,1]:=0.0;" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 42 "for k from 2 to 200 while sol[1,k-1] " 0 "" {MPLTEXT 1 0 2 " " } }{PARA 0 "> " 0 "" {MPLTEXT 1 0 54 " t:=sol[1,k-1]; ylast:=sol[2,k- 1]; h:=sol[3,k-1]; " }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 3 " " }} {PARA 0 "> " 0 "" {MPLTEXT 1 0 25 " #Euler en un paso (h)" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 41 " m1:=deriv(t,ylast);y1:= ylast + h*m1; \+ " }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 1 " " }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 30 " #Euler en dos pasos (h/2) " }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 11 " h2:=h/2;" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 62 " ytemp:=yl ast+h2*m1; m2:=deriv(t+h2,ytemp); y2:=ytemp+h2*m2;" }}{PARA 0 "> " 0 " " {MPLTEXT 1 0 0 "" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 37 " err:=(y2-y 1); abserr:=abs(err); " }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 9 " Uh:=U ;" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 3 " " }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 54 " F:=sqrt(Uh/a bserr); F:=min(2.0,F); F:=max(0.125,F);" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 16 " hnext:=F*h; " }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }} {PARA 0 "> " 0 "" {MPLTEXT 1 0 3 " " }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 45 " sol[1,k]:=sol[1,k-1]+h; sol[2,k]:=y2+err; " }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 34 " sol[3,k]:=hnext; sol[4,k]:=err;" }}{PARA 0 "> \+ " 0 "" {MPLTEXT 1 0 0 "" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 86 " #prin tf(\"K %d -> %5.2f %5.2f (%3.1f) %5.2e %5.2e %5.2e\\n\",k,t+h,hnext,F, y1,y2,err);" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 3 "od;" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 58 "Np:=k-1; \+ printf(\"Numero de pasos en adaptativo %d\\n\",Np);" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 42 "RETURN(linal g[submatrix](sol,1..4,1..Np));" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 4 "end:" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 83 "Vamos a ver como fun ciona nuestro nuevo m\351todo, comparandolo con el Euler original." }} {PARA 0 "" 0 "" {TEXT -1 0 "" }}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 27 "t0:=0.0: tf:=8: y0:= ok(t0):" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 26 "# Soluci\363n con adaptativo:" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 8 "h0:=0.1:" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 8 "E:=1e-3: " }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 34 "adap:=euler_adapta(t0,tf,y0,h0, E):" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 25 "#soluci\363n con Euler fijo:" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 28 "Nfijo:=35: h:=(tf-t0)/Nfijo:" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 38 "eul:=diffnum(ecu,t0,tf,y0,h,foreuler):" }}{PARA 0 "> \+ " 0 "" {MPLTEXT 1 0 0 "" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 39 "#Pintamo s ambos resultados lado a lado." }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 71 " estilo:= style=[point,line]: colores:= color=[red,green]: T:= t=t0..tf :" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 24 "# Adaptativo frente a ok" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 67 "p1:=plot([graf(adap,2),ok(t)],T,estilo,colores,title=`Adaptativo`):" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 16 "# Fijo versus ok" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 29 "colores:= color=[blue,green]:" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 65 "p2:=plot([graf(eul,2),ok(t)],T,estilo,colores,titl e=`Paso Fijo`):" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 29 "display(array(1..2,[p2,p1]));" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 40 "#Comparamos \+ los errores de ambos m\351todos" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 33 " err_eul:=plot_error(eul,ok,blue):" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 34 "err_adap:=plot_error(adap,ok,red):" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 35 "errores:=display(err_eul,err_adap):" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 41 "#Comparamos los pasos usado en adaptativo" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 27 "color es:= color=[red,blue]:" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 69 "pasos:=pl ot([graf(adap,3),h],T,colores,title=`Evolucion del Paso h`):" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 37 "dis play(array(1..2,[errores,pasos]));" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}{PARA 0 "" 0 "" {TEXT -1 1 " " }}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}{PARA 0 "" 0 "" {TEXT -1 232 "Gr\341ficas superio res: Se observa como el m\351todo adaptativo concentra su trabajo cuan do hay zonas complicadas (cambios en el comportamiento de la soluci \363n), como puede verse en la gr\341fica de la derecha. Al usar un pa so fijo nos vemos " }}{PARA 0 "" 0 "" {TEXT -1 124 "forzados a usar un paso innecesariamente peque\361o en el intervalo (2..8) y a la vez d emasiado grande en el intervalo [0..1]." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 191 "Gr\341ficas inferiores: En la gr \341fica de la izquierda se compara el error y en la de la derecha la \+ evoluci\363n del tama\361o del paso: Euler fijo de color azul y el m \351todo adaptativo de color rojo. " }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 144 "Nota 1.- A partir de estas gr\341ficas se puede ver como la estrategia fundamental de un m\351todo adaptativ o es mantener constante el error por paso." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 199 "Nota 2.- Las diferencias entre ambos m\351todos se hacen m\341s notables si alargamos el tiempo de o bservaci\363n, p.e. poner tf=20, 30, 50 y observar el ahorro computaci onal que esta estrategia proporciona. " }}}{SECT 1 {PARA 3 "" 0 "" {TEXT -1 0 "" }{TEXT 256 67 " 6.- Otra estrategia adaptativa: comparar m\351todos de orden p y p+1." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }} {PARA 0 "" 0 "" {TEXT -1 320 " Otra estrategia t\355pica adaptativa es comparar dos m\351todos, uno de orden p y otro de orden t\355picament e p+1. La idea es la misma: ejecutamos los dos m\351todos en cada paso y con las dos aproximaciones de la soluci\363n obtenemos una estimaci \363n del error cometido, y en base a dicha estimaci\363n modificaremo s el tama\361o del paso." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 23 "# Ecuacion diferencial " }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 40 "ecu:=diff(y(t),t)=1/(1+t*t)-2*y(t)*y(t):" } }{PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 20 "# Soluci\363n verdadera" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 35 "ok:= proc(t) RETURN(t/(1+t*t)); end:" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 0 " " }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 35 "# Condiciones iniciales e inter valo" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 21 "t0:=0:tf:=10:y0:=0.0:" }} {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 13 "# y' = f(x,y)" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 50 "deriv := proc \+ (x,y) RETURN(1/(1+x*x)-2*y*y); end:" }{TEXT -1 1 " " }{MPLTEXT 1 0 0 "" }}}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 0 "" }} {PARA 0 "" 0 "" {TEXT -1 0 "" }}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 33 "rk12:= proc(t0,tf,yini,h0,umbral)" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 62 "local t,k,h,ylast,ynext,y1,m1,m2,err,sol,Uh,Np,abserr,F,hnext;" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 25 "sol:=array(1..4,1..1000);" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 56 " sol[1,1]:=t0; sol[2,1]:=y0; sol[3,1]:=h0; sol[4,1]:=0.0;" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 43 "for k f rom 2 to 1000 while sol[1,k-1] " 0 "" {MPLTEXT 1 0 0 "" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 52 " t:=sol[1,k-1]; ylast:=s ol[2,k-1]; h:=sol[3,k-1];" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 4 " " }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 37 " m1:=deriv(t,ylast); y1:=ylast+h *m1;" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 47 " m2:=deriv(t+h,y1); ynext :=ylast+h*(m1+m2)/2;" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 56 " #m2:=deri v(t+0.5*h,ylast+0.5*h*m1); ynext:=ylast+h*m2;" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 18 " err:=(ynext-y1);" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 20 " abserr:=abs(err); " }} {PARA 0 "> " 0 "" {MPLTEXT 1 0 16 " #Uh:=umbral*h;" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 13 " Uh:=umbral;" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 55 " #if (abserr > 2*Uh) then F:= min(2.0,sqrt(Uh/abserr));" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 60 " #eli f (abserr< 0.5*Uh) then F:=max(0.125,sqrt(Uh/abserr));" }}{PARA 0 "> \+ " 0 "" {MPLTEXT 1 0 15 " #else F:=1.0;" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 6 " #fi;" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 1 " " }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 46 " F:= max (0.125, min(2.0, sqrt(Uh/abserr) ));" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 15 " hnext:=F*h; " }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 6 " " }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 65 " so l[1,k]:=t+h; sol[2,k]:=ynext; sol[3,k]:=hnext: sol[4,k]:=err;" }} {PARA 0 "> " 0 "" {MPLTEXT 1 0 68 " #printf(\"%3d %5.2f %5.2f %5.2e % 5.2e %5.2e\\n\",k,t+h,y1,ynext,err);" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 2 " " }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 3 "od;" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 10 " Np:=k-1; " }} {PARA 0 "> " 0 "" {MPLTEXT 1 0 43 " printf(\"Numero de pasos en adap % d\\n\",Np);" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 23 " #for k from 1 to Np do" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 53 " # printf(\"%3d %4.2f %5.2 e \",k,sol[1,k],sol[2,k]); " }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 47 " # p rintf(\"%4.2f %5.2e\\n\",sol[3,k],sol[4,k]); " }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 5 " #od;" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 43 " RETURN(linalg[submatrix](sol,1..4,1..Np) );" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 4 "end:" }}}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 35 "t0:=0: tf:=10: y0:=ok(t0): h0:=0.1:" }}{PARA 0 "> \+ " 0 "" {MPLTEXT 1 0 0 "" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 30 "Nfijo:=2 5: h:=(tf-t0)/Nfijo: " }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 33 "rk2:=diff num(ecu,t0,tf,y0,h,rk2):" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 17 "h0:=0.1 : E:=5e-3:" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 26 "adap:=rk12(t0,tf,y0,h 0,E):" }}}{PARA 0 "" 0 "" {TEXT -1 3 " " }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 39 "#Pintamos ambos resultados lado a lado." }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 71 "estilo:= style=[point,line]: colores:= color= [red,green]: T:= t=t0..tf:" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }} {PARA 0 "> " 0 "" {MPLTEXT 1 0 24 "# Adaptativo frente a ok" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 67 "p1:=plot([graf(adap,2),ok(t)],T,estilo,co lores,title=`Adaptativo`):" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 16 "# Fij o versus ok" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 29 "colores:= color=[blu e,green]:" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 65 "p2:=plot([graf(rk2,2), ok(t)],T,estilo,colores,title=`Paso Fijo`):" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 29 "display(array(1 ..2,[p2,p1]));" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 39 "#Comparamos los errors de ambos m\351todos" }} {PARA 0 "> " 0 "" {MPLTEXT 1 0 33 "err_rk2:=plot_error(rk2,ok,blue):" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 34 "err_adap:=plot_error(adap,ok,red) :" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 35 "errores:=display(err_rk2,err_a dap):" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 41 "#Comparamos los pasos usado en adaprativo" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 27 "colores:= color=[red,blue]:" }}{PARA 0 "> \+ " 0 "" {MPLTEXT 1 0 69 "pasos:=plot([graf(adap,3),h],T,colores,title=` Evoluci\363n del Paso h`):" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }} {PARA 0 "> " 0 "" {MPLTEXT 1 0 37 "display(array(1..2,[errores,pasos]) );" }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}}{PARA 0 "" 0 "" {TEXT -1 3 " \+ " }}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}{PARA 0 "" 0 "" {TEXT -1 1 " " }{MPLTEXT 1 0 1 " " }}}{SECT 1 {PARA 3 "" 0 "" {TEXT -1 0 "" }{TEXT 256 75 " 7.- Ventajas de los m\351todos impl\355citos: \+ estrategia de Predictor-Corrector." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }} {PARA 0 "" 0 "" {TEXT -1 115 " El interes de usar m\351todos impl\355c itos nos lleva a la cuesti\363n de tener que resolver ecuaciones itera tivamente para" }}{PARA 0 "" 0 "" {TEXT -1 115 "hallar el siguiente y_ n. Una estrategia muy com\372n consiste en usar un m\351todo expl\355c ito (predictor) para hallar una " }}{PARA 0 "" 0 "" {TEXT -1 95 "buena hipotesis inicial que luego se refina iterativamente con el m\351todo expl\355cito (corrector)." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 59 "Vamos a programar el siguiente m\351todo predicto r-corrector (" }{XPPEDIT 18 0 "PC^m;" "6#)%#PCG%\"mG" }{TEXT -1 159 " \+ ): Como Predictor un Runge-Kutta de orden 2 (expl\355cito de un paso) \+ y como Corrector un Crank-Nicolson (m\351todo del trapecio) de orden 2 (impl\355cito de un paso)." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 3 " " }}{PARA 0 "" 0 "" {TEXT -1 73 "Estudiaremos la actuaci\363n del m\351todo PC sobre nuesta ecuaci\363n de siempre: " }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 28 "# Nuestra conocida ecuaci\363n." }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 53 "deriv := proc (x,y) RETURN(1/(1+x*x) -2*y*y ); end:" }} {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 39 "ok:=proc(t) RETURN (t/(1+t*t)); end:" }}}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 91 "Aqui definimos el m\351todo de \+ paso fijo h. Adem\341s de los par\341metros habituales hay uno nuevo: " }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 125 " - \+ DELTA, criterio para detener la iteraci\363n. Si iteraciones sucesivas se diferencian menos que DELTA se detiene el proceso." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}} {EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 67 "# Predictor corrector usando RK2 (explicito) y Trapecio (implicito)" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 38 "pred_corr:= proc(t0,tf,y ini,h,m,delta)" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 51 "local tl,t,k,N,so l,y,pred,corr,next,j,dif,yp,m1,m2;" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 22 "N:=round((tf-t0)/h)+1;" }} {PARA 0 "> " 0 "" {MPLTEXT 1 0 27 "tl:=[seq(t0+k*h,k=0..N-1)];" }} {PARA 0 "> " 0 "" {MPLTEXT 1 0 13 "t:=array(tl);" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 29 "yp:=array(1..N) ; yp[1]:=yini;" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 27 "y:=array(1..N); y [1]:=yini;" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 13 "printf(\"\\n\");" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 22 "for k from 1 to N-1 do" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }} {PARA 0 "> " 0 "" {MPLTEXT 1 0 19 "#Prediccion con RK2" }}{PARA 0 "> \+ " 0 "" {MPLTEXT 1 0 63 "#pred:=y[k] + h * deriv(t[k]+h/2, y[k]+h/2*de riv(t[k],y[k]) );" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}{PARA 0 "> \+ " 0 "" {MPLTEXT 1 0 31 "#Prediccion con RK2 (version 2)" }}{PARA 0 "> \+ " 0 "" {MPLTEXT 1 0 24 " m1:=deriv(t[k],y[k]):" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 31 " m2:=deriv(t[k]+h,y[k]+h*m1):" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 31 " pred:=y[k] + h*(m1+m2)/2.0;" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 17 " yp[k+1]:=pre d;" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 75 " printf(\"t=%4.1f OK%7.4f \+ PRED %7.4f -> CORR: \",t[k+1],ok(t[k+1]),pred);" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 28 "#Correccion con trapecio " }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 25 " corr:=pred; dif :=0.1; " }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 41 " for j from 1 to m whi le dif > delta do" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 67 " next:= y[ k] + 0.5*h*(deriv(t[k],y[k]) + deriv(t[k+1],corr) ); " }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 25 " dif:=abs(corr-next);" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 16 " corr:=next;" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 28 " printf(\"%7.4f\",corr); " }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 6 " od;" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 16 " printf(\"\\n\");" } }{PARA 0 "> " 0 "" {MPLTEXT 1 0 16 " y[k+1]:=corr;" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 3 "od;" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }} {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 66 "sol:=array(1..3,1..N,[tl,[seq(yp[k],k=1..N)],[seq(y[k],k=1..N)]]); " }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 22 "#printf(\"N = %d\\n\",N);" }} {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 12 "RETURN(sol);" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 4 "end:" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}} {PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 27 "Aqu\355 e jecutamos el m\351todo (" }{XPPEDIT 18 0 "PC^m;" "6#)%#PCG%\"mG" } {TEXT -1 3 "). " }}{PARA 0 "" 0 "" {TEXT -1 1 " " }}{PARA 0 "" 0 "" {TEXT -1 55 "El listado siguiente detalla las operaciones por pasos:" }}{PARA 0 "" 0 "" {TEXT -1 87 " - Tras PRED mostramos el valor p redicho por el m\351todo expl\355cito rk2 (predictor)." }}{PARA 0 "" 0 "" {TEXT -1 97 " - Tras CORR mostramos las sucesivas iteracion es del corrector. El corrector se para tras " }{XPPEDIT 18 0 "m;" "6# %\"mG" }{TEXT -1 74 " iteraciones o cuando la diferencia entre iterac iones es menor que DELTA:" }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 29 "#Inicializamos las variables." }} {PARA 0 "> " 0 "" {MPLTEXT 1 0 37 "t0:=0.0: y0:=ok(t0): tf:=8.0: h:=0. 5:" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 20 "m:=10: maxdif:=1e-3:" }} {PARA 0 "> " 0 "" {MPLTEXT 1 0 34 "a:=pred_corr(t0,tf,y0,h,m,maxdif): " }}}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 1 " " }{MPLTEXT 1 0 0 "" }{TEXT -1 0 "" }{MPLTEXT 1 0 0 "" } {TEXT -1 147 "Observese que en la zona donde la soluci\363n var\355a m as bruscamente (alrededor del 1.0) el m\351todo corrector ha de ejecut arse un mayor n\372mero de veces." }}{PARA 0 "" 0 "" {TEXT -1 1 " " } {MPLTEXT 1 0 0 "" }{TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 30 "Aqu \355 pintamos los resultados. " }}{PARA 0 "" 0 "" {TEXT -1 0 "" }} {PARA 0 "" 0 "" {TEXT -1 177 "En la gr\341fica de la izquierda mostram os los resultados de ejecutar un RK de orden dos (el predictor). Clara mente para el paso escogido hay problemas en la zona alrededor del 1. \+ " }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 225 "En l a gr\341fica de la derecha mostramos los resultados del predictor-corr ector: en verde la soluci\363n real, en rojo las predicciones efectuad as y en azul la ultima iteraci\363n del corrector (valor definitivo de nuestra soluci\363n). " }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 38 "rk2_solo:=diffnum(ecu,t0,tf,y0,h,rk 2):" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 81 "g1:=plot([ok(t),graf(rk2_sol o,2)],t=t0..tf,style=[line,point],color=[green,red]):" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 95 "g2:=plot([ok(t),graf(a,2),graf(a,3)],t=t0..tf,st yle=[line,point,point],color=[green,red,blue]):" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 29 "display(array(1..2,[g1,g2]));" }}}{EXCHG {PARA 0 "> \+ " 0 "" {MPLTEXT 1 0 0 "" }}}{PARA 0 "" 0 "" {TEXT -1 1 " " }}{PARA 0 " " 0 "" {TEXT -1 111 "Claramente el m\351todo impl\355cito muestra una \+ clara superioridad sobre el expl\355cito a igualdad de orden (ambos so n" }}{PARA 0 "" 0 "" {TEXT -1 115 "de orden 2). Adem\341s, la estrateg ia de predecir el valor permite reducir enormemente el n\372mero de it eraciones para " }}{PARA 0 "" 0 "" {TEXT -1 44 "estimar el valor fina l del m\351todo impl\355cito." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }} {PARA 0 "" 0 "" {TEXT 265 10 "Ejercicio:" }{TEXT -1 108 " 1.- Probar \+ a disminuir el tama\361o del paso h \277tiene alguna ventaja el m\351t odo impl\355cito sobre el expl\355cito?." }}{PARA 0 "" 0 "" {TEXT -1 86 " 2.- Probar a disminuir el valor de MAXDIF aumen tanto el valor de m." }}}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 8 " " }{TEXT 266 2 " " }{TEXT 256 31 "LABORATORIO DE CALCUL O NUMERICO" }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 3 "" 0 "" {TEXT -1 53 " Resoluci\363n num\351rica de ecuaciones diferenciales" }} {PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 92 "En esta pr\341ctica exploraremos algunos de los conceptos b\341sicos en la resoluci\363n num\351rica de " }}{PARA 0 " " 0 "" {TEXT -1 95 "ecuaciones diferenciales, como pueden ser la aprox imaci\363n de la soluci\363n num\351rica a la soluci\363n" }}{PARA 0 " " 0 "" {TEXT -1 60 "exacta y la velocidad de esa aproximaci\363n: orde n del m\351todo." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{EXCHG {PARA 0 "> \+ " 0 "" {MPLTEXT 1 0 60 "restart:with(plots): with(DEtools): with(linal g): Digits=20:" }}{PARA 7 "" 1 "" {TEXT -1 35 "Warning, new definition for adjoint" }}{PARA 7 "" 1 "" {TEXT -1 32 "Warning, new definition f or norm" }}{PARA 7 "" 1 "" {TEXT -1 33 "Warning, new definition for tr ace" }}}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{SECT 1 {PARA 4 "" 0 "" {TEXT -1 53 " 0.- Rutinas auxiliares a usar durante esta pr\341ctica." }} {PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 88 " Maple, p or sus caracteristicas de software matem\341tico eminentemente simb \363lico, tiende a" }}{PARA 0 "" 0 "" {TEXT -1 91 "enmascarar el carac ter num\351rico de algunas de sus funciones. Este es el caso de la fun ci\363n " }}{PARA 0 "" 0 "" {TEXT -1 94 "dsolve( ) que es la que resue lve ecuaciones diferenciales. Dada la naturaleza de nuestro curso" }} {PARA 0 "" 0 "" {TEXT -1 91 "nos interesa dejar bien claro que la sali da de un m\351todo num\351rico de resoluci\363n de ecuacio-" }}{PARA 0 "" 0 "" {TEXT -1 89 "nes diferenciales debe ser un mont\363n de nume ritos, los valores y_n que el m\351todo cree que" }}{PARA 0 "" 0 "" {TEXT -1 88 "son una buena aproximaci\363n a la verdadera ( y desconoc ida) soluci\363n en los puntos t_n. " }}{PARA 0 "" 0 "" {TEXT -1 0 " " }}{PARA 0 "" 0 "" {TEXT -1 92 "Por ello hemos encapsulado la rutina \+ dsolve( ) de Maple dentro de otra, diffnum( ), de forma" }}{PARA 0 "" 0 "" {TEXT -1 89 "que los par\341metros que recibe y devuelve son muy \+ similares a la notaci\363n usada en clase. " }}{PARA 0 "" 0 "" {TEXT -1 92 "El nuevo procedimiento recibe una ecuacion diferencial ecu, un intervalo [t0 tf], un valor " }}{PARA 0 "" 0 "" {TEXT -1 89 "inicial \+ y(t0) y un paso h, as\355 como el nombre de un m\351todo a elegir. Los posibles m\351todos" }}{PARA 0 "" 0 "" {TEXT -1 59 "son los mismos qu e admite dsolve( ) en su opci\363n classical:" }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 5 "" 0 "" {TEXT -1 81 " M\351todos: foreule r, heunform, impoly, rk2, rk3, rk4,adambash,abmoulton, etc" }{TEXT 19 0 "" }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 94 "El procedimiento devuelve una matriz cuya primera fila son los tiempo s tn=t0+n*h, y la segunda" }}{PARA 0 "" 0 "" {TEXT -1 51 "los 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,method=classic al[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 " sol:=arra y(1..2,1..N+1,[ttl,[yy]]);" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 20 " RET URN(eval(sol));" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 4 "end:" }}}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 93 "Al devolver los r esultados de una forma no standard, podemos tener problemas para plote arlos." }}{PARA 0 "" 0 "" {TEXT -1 96 "La siguiente rutina recibe una \+ matriz generada por la rutina anterior y devuelve una estructura " }} {PARA 0 "" 0 "" {TEXT -1 47 "que se puede mandar directamente a un pl ot( )." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 58 "# Rutina que recibe una matriz de dos filas y prepara para" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 53 "# plotear la fila fila 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 88 "Finalmente en los casos en que conozcamos la soluci\363n \+ exacta querremos compararla con la" }}{PARA 0 "" 0 "" {TEXT -1 85 "obt enida num\351ricamente. La siguiente rutina recibe una matriz generada por diffnum( )" }}{PARA 0 "" 0 "" {TEXT -1 98 "(formada por tn y yn) \+ y una funci\363n f( ), la soluci\363n exacta. La rutina calcula en lo s diferentes" }}{PARA 0 "" 0 "" {TEXT -1 64 "tn la discrepancia entre f(tn) e yn, devolviendonos un gr\341fico." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 64 "# Recibimos una so luci\363n numerica (matriz a) y solucion ok (fun)" }}{PARA 0 "> " 0 " " {MPLTEXT 1 0 0 "" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 27 "plot_error:=p roc(a,fun,col)" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 21 " local t,y,g,k,d ata;" }}{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 58 " data:=[seq([t[k ],abs( y[k]-fun(t[k]))],k=1..coldim(a))];" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 3 " " }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 40 " g:=plot(da ta,color=col,title=`Error`):" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 12 " R ETURN(g);" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 4 "end:" }}}{PARA 0 "" 0 " " {TEXT -1 0 "" }}}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{SECT 1 {PARA 3 "" 0 " " {TEXT -1 1 " " }{TEXT 256 63 "1.- Aproximaci\363n de la soluci\363n \+ num\351rica a la soluci\363n exacta. " }}{PARA 3 "" 0 "" {TEXT 260 0 " " }{TEXT 258 0 "" }{TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 87 "Ahora \+ nos planteamos la cuesti\363n principal que nos surge al resolver num \351ricamente una" }}{PARA 0 "" 0 "" {TEXT 267 21 "ecuaci\363n diferen cial:" }}{PARA 0 "" 0 "" {TEXT 262 103 "Cual es la aproximaci\363n de \+ la soluci\363n num\351rica a la soluci\363n exacta. Veamos un ejemplo \+ ilustrativo. " }}{PARA 0 "" 0 "" {TEXT -1 1 " " }}{PARA 0 "" 0 "" {TEXT 263 49 "El ejemplo del cazador-presa con diversos m\351todos" }} {PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 91 "En este a partado resolveremos de nuevo el problema cazador-presa pero con diver sos m\351todos:" }}{PARA 0 "" 0 "" {TEXT -1 87 "Euler, Runge-Kutta de \+ diversos ordenes, etc. y diversos tama\361os de paso. Todo ello lo " }}{PARA 0 "" 0 "" {TEXT -1 94 "compararemos con la soluci\363n \"exact a\", las \363rbitas peri\363dicas que obteniamos con el m\351todo por " }}{PARA 0 "" 0 "" {TEXT -1 22 "defecto que usa Maple." }}{PARA 0 "" 0 "" {TEXT -1 1 " " }}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}} {EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 36 "# Constantes que definen el \+ problema" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 38 "a:=0.6: # Crecimie nto de roedores" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 66 "b:=-0.01: # D ecaimiento de roedores al encontrarse con un zorro" }}{PARA 0 "> " 0 " " {MPLTEXT 1 0 60 "c:=-0.5: # DEcaimiento de zorros en ausencia de roedores" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 65 "d:=0.001: # Benefic io de un zorro al encontrarse con un roedor" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 12 "# Ecuaciones" } }{PARA 0 "> " 0 "" {MPLTEXT 1 0 35 "ec1:=D(r)(t) = a*r(t) +b*r(t)*z(t) :" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 36 "ec2:=D(z)(t) = c*z(t) + d*r(t) *z(t):" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 23 "# Condiciones iniciales" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 16 "ini1:= r(0)=600:" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 15 "ini2:= z(0)=40:" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 36 "# Resoluci\363n de la ecua ci\363n \"exacta\"" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 52 "ok:=dsolve(\{ ec1,ec2,ini1,ini2\},\{r(t),z(t)\},numeric):" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}{PARA 0 "" 0 "" {TEXT -1 1 " " }{MPLTEXT 1 0 0 "" }}}{PARA 0 "" 0 "" {TEXT -1 110 "A continuaci\363n resolvemos por u n cierto m\351todo y paso, pintando la soluci\363n \"exacta\" en rojo (linea continua)" }}{PARA 0 "" 0 "" {TEXT -1 65 "y la otra con puntos azules unidos por linea discontinua (verde):" }}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}{PARA 0 "" 0 "" {TEXT -1 1 " " }{MPLTEXT 1 0 0 "" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 62 "# Resoluci\363n de l a ecuaci\363n por diversos m\351todos que conocemos" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 40 "# Especifica ci\363n de m\351todo y paso a usar" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 17 "metodo:=foreuler;" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 7 "h:=0.1;" }} {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 31 "# Tiempo de observaci\363n (meses)" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 9 "Tobs:=36;" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 19 "N:=round(Tobs/h )+1:" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 29 "tt:=array([seq(k*h,k=0..N)] ):" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 98 "sol:=dsolve(\{ec1,ec2,ini1,ini2\},\{z(t),r(t)\},numer ic,method=classical[metodo],stepsize=h,value=tt):" }}{PARA 0 "> " 0 " " {MPLTEXT 1 0 0 "" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 16 "zr:=[z(t),r(t)]:" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 28 "ver:= view=[0..180,0..1200]:" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 48 "g1:=odeplot(ok,zr,0..16,numpoints=50,color=red):" }} {PARA 0 "> " 0 "" {MPLTEXT 1 0 55 "g2:=odeplot(sol,zr,0..Tobs,color=bl ue,style=point,ver):" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 56 "g3:=odeplot (sol,zr,0..Tobs,color=green,linestyle=3,ver):" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 18 "display(g1,g2,g 3);" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#>%'metodoG%)foreulerG" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#>%\"hG$\"\"\"!\"\"" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#>%%TobsG\"#O" }}{PARA 13 "" 1 "" {GLPLOT2D 400 300 300 {PLOTDATA 2 "6&-%'CURVESG6$7U7$$\"#S\"\"!$\"$+'F*7$$\"1fWT,]WbT!#9$\"1 %z\\(eKh\"Q'!#87$$\"19lr)H!3pVF0$\"12l$oOAwu'F37$$\"1b#*3>%*\\XYF0$\"1 exTn$y*yqF37$$\"1Pd&4!>$y)\\F0$\"1GjE7G[`tF37$$\"1YE?25-'R&F0$\"1%z\"[ LjSZvF37$$\"1t#p%\\VgkeF0$\"1.-SI4NQwF37$$\"1g0)Qa3/Q'F0$\"1!Q8!3>B4wF 37$$\"1.G1Dw+@pF0$\"15[%y>OEX(F37$$\"1-)>kD!3buF0$\"1=!*4ok:urF37$$\"1 GT**oz_XzF0$\"1C7uRP5$z'F37$$\"1?kz')QswQ=#))F0$\"1-xozR+d`F37$$\"1A_)3**4f&) )F0$\"1BB6L%z+*[F37$$\"1BH!y7FUw)F0$\"1:\"*)G1O\"oWF37$$\"1wGb\"*Gkk&) F0$\"1Y?Lj;!>5%F37$$\"1a\"3c$)\\*z#)F0$\"1!)Hw<*>_z$F37$$\"1@kC5W0MzF0 $\"1U#z5u(HZNF37$$\"1A5W\"y+$\\vF0$\"1)=\"Q1muaLF37$$\"1i0WTs)\\9(F0$ \"1BY;KF37$$ \"1mTbTK?+ZF0$\"1%GVcyuDM$F37$$\"1(yxAZ5\"oWF0$\"1.+wM<#y\\$F37$$\"1zF #>=[3F%F0$\"1@L*3#o^&o$F37$$\"1q-VC:-4TF0$\"1N%ypZlc!RF37$$\"1xw/SRL$) RF0$\"1'Rn8Yd!eTF37$$\"1)y>r\"*zZ*QF0$\"1X#R;1u>W%F37$$\"1;1QiT\"[%QF0 $\"1'*eFsu#ev%F37$$\"1jn%4s(\\NQF0$\"1CmSr7s'4&F37$$\"1/$o@D7'pQF0$\"1 5s5L%H*faF37$$\"1`eDl!>2&RF0$\"1(*)46Y4$QeF37$$\"1Dp:!ffjrF;F%F0$\"1EjxqZR'f'F37$$\"10*4:<&3@XF0$\"1K(R/P0\\%pF 37$$\"1bVW>DTN[F0$\"1Tt(z\"yAYsF37$$\"1H0;cF9;_F0$\"1*z)R%Q.qZ(F37$$\" 1A\\8XCQgcF0$\"1_r?QI%Qh(F37$$\"1Qv4uZYehF0$\"1v4^xu)oj(F37$$\"1HG3#4c >p'F0$\"1&[e7gVT`(F37$$\"1xx#\\cuIB(F0$\"1!z:u7^aI(F37$$\"1o^X#4!fYxF0 $\"1GyfZ=ZkpF37$$\"12Chw!\\X>)F0$\"1(Q(=S4SPlF37$$\"1qlJ3E#Ga)F0$\"1(* >Sl3NegF37$$\"1^rlxbJn()F0$\"1Nnk#)o#Hc&F37$$\"1]nN#>Mw&))F0$\"1p:*p$* 3B3&F3-%'COLOURG6&%$RGBG$\"*++++\"!\")F*F*-F$6%7falF'7$$\"++++SSFi[l$ \"++++?h!\"(7$$\"++![_3%Fi[l$\"++?&*RiFb\\l7$$\"+9J!f8%Fi[l$\"+1;VfjFb \\l7$$\"+))y7#>%Fi[l$\"+tv(zZ'Fb\\l7$$\"+)f'3aUFi[l$\"+B64&f'Fb\\l7$$ \"+cJ%>K%Fi[l$\"+.dB5nFb\\l7$$\"+r&eeR%Fi[l$\"+ts$G#oFb\\l7$$\"+?z)fZ% Fi[l$\"+K_GKpFb\\l7$$\"+unZiXFi[l$\"+'4Mz.(Fb\\l7$$\"+WqXbYFi[l$\"+Lg5 RrFb\\l7$$\"+2A/bZFi[l$\"+!Q%4NsFb\\l7$$\"+))=Kh[Fi[l$\"+^#o^K(Fb\\l7$ $\"+,eNu\\Fi[l$\"+R$y&3uFb\\l7$$\"+dq;%4&Fi[l$\"+$*Rc%[(Fb\\l7$$\"+(*[ t?_Fi[l$\"+R;O_vFb\\l7$$\"+Kp)RN&Fi[l$\"+eX@6wFb\\l7$$\"+=5z$\\&Fi[l$ \"+**RQgwFb\\l7$$\"+bp%*RcFi[l$\"+^:;*p(Fb\\l7$$\"+N#y@z&Fi[l$\"+*RM(Fi[l$\"+jx]4uF b\\l7$$\"+!)f%4_(Fi[l$\"+xX#*4tFb\\l7$$\"+gTn%p(Fi[l$\"+tXu)>(Fb\\l7$$ \"+\"RgQ'yFi[l$\"+4$\\n2(Fb\\l7$$\"+dItP+lFb\\l7$$\"+1_!ff)Fi[l$\"+iL*)RjFb\\l7$$\"+p636()Fi[l$\"+SdJvh Fb\\l7$$\"+xQY8))Fi[l$\"+=z*y+'Fb\\l7$$\"+'f%H-*)Fi[l$\"+!)y')QeFb\\l7 $$\"+'3tp(*)Fi[l$\"+LnSpcFb\\l7$$\"+[b1P!*Fi[l$\"+?+j+bFb\\l7$$\"+6yI# 3*Fi[l$\"+\"GsNL&Fb\\l7$$\"+sQg7\"*Fi[l$\"+m^;&[Fb\\l7$$\"+<0B:\"*Fi[l$\"+YH #)*p%Fb\\l7$$\"+(ooy3*Fi[l$\"+FET`XFb\\l7$$\"+8MGZ!*Fi[l$\"+C#4GT%Fb\\ l7$$\"+*feT**)Fi[l$\"+?%Q$yUFb\\l7$$\"+77DH*)Fi[l$\"+#=Q-:%Fb\\l7$$\"+ XQP`))Fi[l$\"+!Gn'GSFb\\l7$$\"+G\"ytw)Fi[l$\"+UVr8RFb\\l7$$\"+d$R@n)Fi [l$\"+nqS0QFb\\l7$$\"+$fU&o&)Fi[l$\"+(G@Pq$Fb\\l7$$\"+*RquX)Fi[l$\"+G' *e3OFb\\l7$$\"+zAzR$)Fi[l$\"+1'4*>NFb\\l7$$\"+:eN;#)Fi[l$\"+L5bPMFb\\l 7$$\"+[%zz3)Fi[l$\"+rEOhLFb\\l7$$\"+)zYa&zFi[l$\"+4\"y6H$Fb\\l7$$\"+#Q -&>yFi[l$\"+s3#oA$Fb\\l7$$\"+)f[3o(Fi[l$\"+*y3\"oJFb\\l7$$\"+3Q9SvFi[l $\"+xw&[6$Fb\\l7$$\"+m8+)R(Fi[l$\"+*Q%)o1$Fb\\l7$$\"+Z%*)\\D(Fi[l$\"+. $4S-$Fb\\l7$$\"+I:j6rFi[l$\"+0$eg)HFb\\l7$$\"+ikGFb\\l7$$\"+kOBoiFi[l$\"+1 ,-`GFb\\l7$$\"+[flLhFi[l$\"+YtOXGFb\\l7$$\"+%>)\\,gFi[l$\"+VVcTGFb\\l7 $$\"+<(f>(eFi[l$\"+r<_TGFb\\l7$$\"+XZ@XdFi[l$\"+k+;XGFb\\l7$$\"+jXT@cF i[l$\"+7\"4C&GFb\\l7$$\"+w&*o+bFi[l$\"+=z?jGFb\\l7$$\"+x7:$Q&Fi[l$\"+9 U]xGFb\\l7$$\"+=T*)o_Fi[l$\"+iSD&*GFb\\l7$$\"+$G(*z:&Fi[l$\"+N9U;HFb\\ l7$$\"+dj_]]Fi[l$\"+$yx4%HFb\\l7$$\"+(*[`Y\\Fi[l$\"+#f,*oHFb\\l7$$\"+0 f1Y[Fi[l$\"+Mz<+IFb\\l7$$\"+5K:\\ZFi[l$\"+5!)zMIFb\\l7$$\"+dF#el%Fi[l$ \"+%oeF2$Fb\\l7$$\"+,)*zWFi[l$\"+NZreJFb \\l7$$\"+04\\(R%Fi[l$\"+Fys1KFb\\l7$$\"+C>j=VFi[l$\"+Kf6eKFb\\l7$$\"+j jSVUFi[l$\"+`o*GJ$Fb\\l7$$\"+@d\"=<%Fi[l$\"+))44rLFb\\l7$$\"+R1'Q5%Fi[ l$\"+W2sKMFb\\l7$$\"+$oT&RSFi[l$\"+8*4y\\$Fb\\l7$$\"+)3g)yRFi[l$\"+=IQ mNFb\\l7$$\"+u%==#RFi[l$\"+4YYQOFb\\l7$$\"+G:UoQFi[l$\"+3&ySr$Fb\\l7$$ \"+kmn=QFi[l$\"+1qC$z$Fb\\l7$$\"+pYfsPFi[l$\"+))**)f(QFb\\l7$$\"+F.>IP Fi[l$\"+'*RKiRFb\\l7$$\"+PI[\"p$Fi[l$\"+47E_SFb\\l7$$\"+?u\\cOFi[l$\"+ Y$3e9%Fb\\l7$$\"+CREDOFi[l$\"+sa'HC%Fb\\l7$$\"+C%>yf$Fi[l$\"+/ZsVVFb\\ l7$$\"+Cy?uNFi[l$\"+:)o![WFb\\l7$$\"+c1[aNFi[l$\"+?(pfb%Fb\\l7$$\"+'o( pQNFi[l$\"+ToQnYFb\\l7$$\"+7v#p_$Fi[l$\"+T`E#y%Fb\\l7$$\"+o\"[#>NFi[l$ \"+IU`+\\Fb\\l7$$\"+@xu:NFi[l$\"+BV5A]Fb\\l7$$\"+`[_;NFi[l$\"+kg'o9&Fb \\l7$$\"+N%*o@NFi[l$\"+/sou_Fb\\l7$$\"+lIOJNFi[l$\"+M.T0aFb\\l7$$\"+w& zca$Fi[l$\"+!H])QbFb\\l7$$\"+\"[&ykNFi[l$\"+C9zucFb\\l7$$\"+[.%))e$Fi[ l$\"+pZ)H\"eFb\\l7$$\"+tq,=OFi[l$\"+7^9`fFb\\l7$$\"+E?]_OFi[l$\"+7![\\ 4'Fb\\l7$$\"+C]\\#p$Fi[l$\"+$zE!QiFb\\l7$$\"+-\"4#QPFi[l$\"+s&p>Q'Fb\\ l7$$\"+9+()*y$Fi[l$\"+zjJElFb\\l7$$\"+@arZQFi[l$\"+akbqmFb\\l7$$\"+zO* >\"RFi[l$\"+7e79oFb\\l7$$\"+m@'H)RFi[l$\"+!>0k&pFb\\l7$$\"+X])31%Fi[l$ \"+:&=n4(Fb\\l7$$\"+u..YTFi[l$\"+X?LMsFb\\l7$$\"+fkmQUFi[l$\"+kVXotFb \\l7$$\"+St0RVFi[l$\"+AuB)\\(Fb\\l7$$\"+%HduW%Fi[l$\"+Y)yFi(Fb\\l7$$\" +jU5kXFi[l$\"+Vd7TxFb\\l7$$\"++@@*o%Fi[l$\"+Q-G_yFb\\l7$$\"+q:'H#[Fi[l $\"+xp?bzFb\\l7$$\"+n+\\l\\Fi[l$\"+?G%)[!)Fb\\l7$$\"+a+)o6&Fi[l$\"+**) 3@8)Fb\\l7$$\"+5j9x_Fi[l$\"+v^#R?)Fb\\l7$$\"+H@AYaFi[l$\"+^vAj#)Fb\\l7 $$\"+]Z%Ri&Fi[l$\"+xu)*3$)Fb\\l7$$\"+y./5eFi[l$\"+hQBS$)Fb\\l7$$\"+9$4 T+'Fi[l$\"+Qp2c$)Fb\\l7$$\"+U=h0iFi[l$\"+hNtb$)Fb\\l7$$\"+Tc&QT'Fi[l$ \"+$=`&Q$)Fb\\l7$$\"+Bc)zi'Fi[l$\"+5O//$)Fb\\l7$$\"+er(p%oFi[l$\"+5a*= D)Fb\\l7$$\"+%pL'pqFi[l$\"+TP+#=)Fb\\l7$$\"+.(*e%H(Fi[l$\"+si[%4)Fb\\l 7$$\"+!y>._(Fi[l$\"+')ep*)zFb\\l7$$\"+rW:XxFi[l$\"+fpAoyFb\\l7$$\"+!3. t'zFi[l$\"+WU\"4t(Fb\\l7$$\"+&H$)[=)Fi[l$\"+IP#)yvFb\\l7$$\"+3q&fR)Fi[ l$\"+w_B8uFb\\l7$$\"+/7d)f)Fi[l$\"+YthNsFb\\l7$$\"+`B!3z)Fi[l$\"+xYfZq Fb\\l7$$\"+aB!3(*)Fi[l$\"+R-\"4&oFb\\l7$$\"+`Q%o8*Fi[l$\"+OKQZmFb\\l7$ $\"+(oitG*Fi[l$\"+/_')QkFb\\l7$$\"+<`*4U*Fi[l$\"+_j>FiFb\\l7$$\"+G%4m` *Fi[l$\"+cU;9gFb\\l7$$\"+7jKL'*Fi[l$\"+cnY,eFb\\l7$$\"+?U`5(*Fi[l$\"+Q 0o!f&Fb\\l7$$\"+cC*yw*Fi[l$\"+CkB$Q&Fb\\l7$$\"+ylK0)*Fi[l$\"+l=S!=&Fb \\l7$$\"+qb,B)*Fi[l$\"+c1F$)\\Fb\\l7$$\"+WAP@)*Fi[l$\"+V%fFz%Fb\\l7$$ \"+y$=5!)*Fi[l$\"+k-h4YFb\\l7$$\"+3ivi(*Fi[l$\"+J\")RMWFb\\l7$$\"+wy`2 (*Fi[l$\"+SDanUFb\\l7$$\"+\"HMkj*Fi[l$\"+%y@$4TFb\\l7$$\"+qYg]&*Fi[l$ \"+)**)))fRFb\\l7$$\"+!yn7X*Fi[l$\"+%**)G>QFb\\l7$$\"+@cnR$*Fi[l$\"+/^ Z(o$Fb\\l7$$\"+`+4<#*Fi[l$\"+(RDVc$Fb\\l7$$\"+KEw%3*Fi[l$\"+Ryl\\MFb\\ l7$$\"+Fx\"R%*)Fi[l$\"+$3WKM$Fb\\l7$$\"+T)Qdz)Fi[l$\"+G<#[C$Fb\\l7$$\" +WzNT')Fi[l$\"+()\\5aJFb\\l7$$\"+\\v%=[)Fi[l$\"+%y$zqIFb\\l7$$\"+>_@=$ )Fi[l$\"+k8e%*HFb\\l7$$\"+\"=+9:)Fi[l$\"+B0;DHFb\\l7$$\"+7`#Fb\\l7$$\"+_@v^hFi[l$\"+A->CDFb\\l7 $$\"+rkW**fFi[l$\"+4(f._#Fb\\l7$$\"+zyUFi[l$\"+gm6VGFb\\l7$$\"+a0#f=%Fi[l$\"+\\&p?*GFb\\l 7$$\"+gUo(4%Fi[l$\"+))R`WHFb\\l7$$\"+`xX8SFi[l$\"+A$[0+$Fb\\l7$$\"+01@ LRFi[l$\"+#[b,1$Fb\\l7$$\"+NC\"p&QFi[l$\"+^CSBJFb\\l7$$\"+6P`%y$Fi[l$ \"+,(R.>$Fb\\l7$$\"+2l/;PFi[l$\"+-1-hKFb\\l7$$\"+C_U^OFi[l$\"+'z+bL$Fb \\l7$$\"+!GZ1f$Fi[l$\"+Ev$QT$Fb\\l7$$\"+!y$pLNFi[l$\"+8*)3'\\$Fb\\l7$$ \"+o,b![$Fi[l$\"+rJJ#e$Fb\\l7$$\"+no?JMFi[l$\"+axcsOFb\\l7$$\"+@*fcQ$F i[l$\"+A%3pw$Fb\\l7$$\"+M;\"RM$Fi[l$\"+>#)QlQFb\\l7$$\"+>7(fI$Fi[l$\"+ Yj0oRFb\\l7$$\"+fa&=F$Fi[l$\"+Ep&\\2%Fb\\l7$$\"+(Q*eTKFi[l$\"+Tw7'=%Fb \\l7$$\"+()p?:KFi[l$\"+I#)f,VFb\\l7$$\"+C>v#>$Fi[l$\"+P))Q@WFb\\l7$$\" +0$yU<$Fi[l$\"+!>3ba%Fb\\l7$$\"+q9&)fJFi[l$\"+,;&Rn%Fb\\l7$$\"+D)[&\\J Fi[l$\"+p()p1[Fb\\l7$$\"+32YVJFi[l$\"+r8rV\\Fb\\l7$$\"+(H\"pTJFi[l$\"+ H/$\\3&Fb\\l7$$\"+\\&fV9$Fi[l$\"+QMFI_Fb\\l7$$\"+u,g^JFi[l$\"+U7jz`Fb \\l7$$\"+LYcjJFi[l$\"+\\Y'G`&Fb\\l7$$\"+^@U!=$Fi[l$\"+z2!)*o&Fb\\l7$$ \"+G2O-KFi[l$\"+T\"H-&eFb\\l7$$\"+K\")eHKFi[l$\"+\\u*Q,'Fb\\l7$$\"+WGL iKFi[l$\"+yr]!='Fb\\l7$$\"+<\\%3I$Fi[l$\"+$*)3(\\jFb\\l7$$\"+9nRXLFi[l $\"+$Q(4@lFb\\l7$$\"+cMG'R$Fi[l$\"+[m?%p'Fb\\l7$$\"+;N#QX$Fi[l$\"+9[]o oFb\\l7$$\"+v$e$=NFi[l$\"+o!*QVqFb\\l7$$\"+D@D!f$Fi[l$\"+V2==sFb\\l7$$ \"+1/*)pOFi[l$\"++27#R(Fb\\l7$$\"+9'ywv$Fi[l$\"+9_OkvFb\\l7$$\"+<#RS&Q Fi[l$\"+$f#)Rt(Fb\\l7$$\"+yzSfRFi[l$\"+E3&***yFb\\l7$$\"+i)GU2%Fi[l$\" +#fc61)Fb\\l7$$\"+(QZ*)>%Fi[l$\"+>gR;#)Fb\\l7$$\"+k@+MVFi[l$\"+MwPk$)F b\\l7$$\"+kV\")zWFi[l$\"+%)zs.&)Fb\\l7$$\"+u[xOYFi[l$\"+N/+L')Fb\\l7$$ \"+D*G_![Fi[l$\"+mwo]()Fb\\l7$$\"+a!ea)\\Fi[l$\"+^$Q_&))Fb\\l7$$\"+X$f w<&Fi[l$\"+f%y]%*)Fb\\l7$$\"+o?#>Q&Fi[l$\"+wuj=!*Fb\\l7$$\"+0??)f&Fi[l $\"+%ozV2*Fb\\l7$$\"+5SHEeFi[l$\"+g.%36*Fb\\l7$$\"+cO!e1'Fi[l$\"+NkmE \"*Fb\\l7$$\"+S\">hJ'Fi[l$\"+b2m?\"*Fb\\l7$$\"+$)\\QwlFi[l$\"+,'G=4*Fb \\l7$$\"+&Qza%oFi[l$\"+mYUR!*Fb\\l7$$\"+ot*>7(Fi[l$\"+%>)*H'*)Fb\\l7$$ \"+&*>C/uFi[l$\"+!fLC'))Fb\\l7$$\"+DfA!p(Fi[l$\"+wN)zt)Fb\\l7$$\"+2`ox zFi[l$\"+7>H!f)Fb\\l7$$\"+(\\2TE)Fi[l$\"+rHS?%)Fb\\l7$$\"+vKxY&)Fi[l$ \"+'*fvH#)Fb\\l7$$\"+(>8G#))Fi[l$\"+qF;?!)Fb\\l7$$\"+;lF*3*Fi[l$\"+d&o Pz(Fb\\l7$$\"+l)4KM*Fi[l$\"+&\\(*Hb(Fb\\l7$$\"+z/+RqFb\\l7$$\"+cM?+5Fb\\l$\"+wdOrnFb\\l7$$\"+ =2#z,\"Fb\\l$\"+yLP+lFb\\l7$$\"+HL>L5Fb\\l$\"+2$4(GiFb\\l7$$\"+s#))e/ \"Fb\\l$\"+vd))efFb\\l7$$\"+Yr\"f0\"Fb\\l$\"+Zg=$p&Fb\\l7$$\"+;mBj5Fb \\l$\"+MRiLaFb\\l7$$\"+kq%y1\"Fb\\l$\"+y&=>=&Fb\\l7$$\"+w'*yp5Fb\\l$\" +RS[R\\Fb\\l7$$\"+$GU\"p5Fb\\l$\"+y@V2ZFb\\l7$$\"+V5Fb\\l$\"+%)3r&*QFb\\l7$$\"+^ynJ5Fb\\l$\"+YQ0BPFb\\l7$$\"+!)Q] =5Fb\\l$\"+eyLiNFb\\l7$$\"+N7'Q+\"Fb\\l$\"+PED8MFb\\l7$$\"+L\\Kz)*Fi[l $\"+%f/aF$Fb\\l7$$\"+\\l%*3(*Fi[l$\"+25M[JFb\\l7$$\"+n*p\"H&*Fi[l$\"+A 2dJIFb\\l7$$\"+&*\\fT$*Fi[l$\"+R9eCHFb\\l7$$\"+*z_L)Fi[l$\"+e&y%=DFb\\l7$$\"+!yy$G \")Fi[l$\"+olmfCFb\\l7$$\"+c3*=#zFi[l$\"+ZbJ2CFb\\l7$$\"+B7];xFi[l$\"+ q&\\5O#Fb\\l7$$\"+\"emG^(Fi[l$\"+D@_?BFb\\l7$$\"+$)4c6tFi[l$\"+@dT&G#F b\\l7$$\"+)\\#38rFi[l$\"++6WbAFb\\l7$$\"+f(fy\"pFi[l$\"+!=O.B#Fb\\l7$$ \"++$eis'Fi[l$\"+A[')4AFb\\l7$$\"+,meQlFi[l$\"+&\\:Q>#Fb\\l7$$\"+(z,^N 'Fi[l$\"+**)**>=#Fb\\l7$$\"+A],whFi[l$\"+y:Du@Fb\\l7$$\"+vj\\,gFi[l$\" +pXUq@Fb\\l7$$\"+3&z;$eFi[l$\"+\"4#Rq@Fb\\l7$$\"+VemmcFi[l$\"+1`/u@Fb \\l7$$\"+M%Gl]&Fi[l$\"+U@H\"=#Fb\\l7$$\"+%[:8N&Fi[l$\"+*>c?>#Fb\\l7$$ \"+`N0,_Fi[l$\"+GdF1AFb\\l7$$\"+d/vb]Fi[l$\"+!p-RA#Fb\\l7$$\"+0zR:\\Fi [l$\"+\")=!\\C#Fb\\l7$$\"+qQ(*zZFi[l$\"+M,DpAFb\\l7$$\"+SZW\\YFi[l$\"+ zb$pH#Fb\\l7$$\"+aswBXFi[l$\"+jp&zK#Fb\\l7$$\"+Y-*GS%Fi[l$\"+EIKiBFb\\ l7$$\"+Civ'G%Fi[l$\"+=>0+CFb\\l7$$\"+)y-`<%Fi[l$\"+e1XNoQFi[l$\"+! ye_e#Fb\\l7$$\"+oP%\\x$Fi[l$\"+KtOSEFb\\l7$$\"+i*oeo$Fi[l$\"+!*p6*p#Fb \\l7$$\"+H91,OFi[l$\"+%4y:w#Fb\\l7$$\"+\"[a/_$Fi[l$\"+cm#y#GFb\\l7$$\" +2T)RW$Fi[l$\"+0R%z*GFb\\l7$$\"+H'*erLFi[l$\"+7e,sHFb\\l7$$\"+CV@.LFi[ l$\"+%eK,0$Fb\\l7$$\"+Cg!)QKFi[l$\"+B\")QKJFb\\l7$$\"+)p<$yJFi[l$\"+N% z)=KFb\\l7$$\"+G!3<7$Fi[l$\"+')fq4LFb\\l7$$\"+%)>%*oIFi[l$\"+))*o\\S$F b\\l7$$\"+39**>IFi[l$\"+01x/NFb\\l7$$\"+4c$[(HFi[l$\"+qI@4OFb\\l7$$\"+ %)>YLHFi[l$\"+*p(R=PFb\\l7$$\"+om'e*GFi[l$\"+yPUKQFb\\l7$$\"+<_0iGFi[l $\"+AtQ^RFb\\l7$$\"+VL/KGFi[l$\"+u'z`2%Fb\\l7$$\"+&pde!GFi[l$\"+Of[/UF b\\l7$$\"+-n`$y#Fi[l$\"++KyQVFb\\l7$$\"+!\\J^w#Fi[l$\"+p&Q$yWFb\\l7$$ \"+moq]FFi[l$\"+Lp?BYFb\\l7$$\"+$RU.u#Fi[l$\"+z%GMx%Fb\\l7$$\"+eN8MFFi [l$\"+/f-H\\Fb\\l7$$\"+JI>KFFi[l$\"+28+!4&Fb\\l7$$\"+S?lMFFi[l$\"+CFLc _Fb\\l7$$\"+[=mTFFi[l$\"+x-(zU&Fb\\l7$$\"+Y`R`FFi[l$\"+.>$[g&Fb\\l7$$ \"+g([+x#Fi[l$\"+O')z'y&Fb\\l7$$\"+mM%=z#Fi[l$\"+5%4P(fFb\\l7$$\"+4z-> GFi[l$\"+f`NlhFb\\l7$$\"+-'z=&GFi[l$\"+#ft9O'Fb\\l7$$\"+!>22*GFi[l$\"+ S/uhlFb\\l7$$\"+LD&e$HFi[l$\"+kTwlnFb\\l7$$\"+lFp()HFi[l$\"+br2tpFb\\l 7$$\"+[AkYIFi[l$\"+jw7$=(Fb\\l7$$\"+QV:8JFi[l$\"+@6F&R(Fb\\l7$$\"+DGs( =$Fi[l$\"+&=h(3wFb\\l-Fd[l6&Ff[lF*F*Fg[l-%&STYLEG6#%&POINTG-F$6%F\\\\l -Fd[l6&Ff[lF*Fg[lF*-%*LINESTYLEG6#\"\"$-%%VIEWG6$;F*$\"$!=F*;F*$\"%+7F *" 1 2 0 1 0 2 9 1 4 2 1.000000 45.000000 45.000000 0 }}}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}{PARA 0 "" 0 "" {TEXT -1 1 " " } {MPLTEXT 1 0 0 "" }}}{PARA 256 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 " " {TEXT 264 10 "Ejercicio:" }{TEXT -1 1 " " }}{PARA 0 "" 0 "" {TEXT -1 91 "a) Hacer pruebas con Euler con h = 0.2, 0.1, 0.05 y compro bad que tiende a irse hacia " }}{PARA 0 "" 0 "" {TEXT -1 93 "afuera (l \363gico pues Euler sigue la tangente y \351sta siempre nos echa hacia afuera). \277Qu\351 valor" }}{PARA 0 "" 0 "" {TEXT -1 66 "de h se pre cisa para que la soluci\363n se confunda con la verdadera?" }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 84 "b) Probad con rk2 . Ver como con h=0.2 funciona mucho mejor que Euler. Sin embargo no" } }{PARA 0 "" 0 "" {TEXT -1 89 "est\341 libre de problemas. Si ponemos h =0.95 y Tobs = 400 el sistema parece estabilizarse" }}{PARA 0 "" 0 " " {TEXT -1 88 "en una \363rbita m\341s grande. Si s\363lo hubiesemos v isto esta soluci\363n podr\355amos pensar que se" }}{PARA 0 "" 0 "" {TEXT -1 48 " alcanza un equilibrio distinto del \"verdadero\"." }} {PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 66 "c) Probad con rk3 y paso h=1.0, \277hacia d\363nde va ahora la soluci\363n?" }} {PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 103 "d) Final mente mostrad como rk4 es mucho m\341s preciso, manteniendose sobre la soluci\363n incluso para h=1.0" }{MPLTEXT 1 0 0 "" }}}{PARA 0 "" 0 " " {TEXT -1 0 "" }}{SECT 1 {PARA 3 "" 0 "" {TEXT 268 21 "2.- Orden del \+ m\351todo." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 1 " " }}{PARA 0 "" 0 "" {TEXT -1 94 "Dentro de los m\351todos conve rgentes (esto es, aquellos que se aproximan a la verdadera soluci\363n " }}{PARA 0 "" 0 "" {TEXT -1 20 "al disminuir h), el " }{TEXT 256 6 "o rden " }{TEXT -1 64 "de un m\351todo estudia la \"velocidad\" a la que esa aproximaci\363n se" }}{PARA 0 "" 0 "" {TEXT -1 40 "realiza. Veamo s gr\341ficamente este hecho." }}{PARA 0 "" 0 "" {TEXT -1 1 " " }} {PARA 0 "" 0 "" {TEXT -1 95 "Si probamos varios m\351todos y pintamos la evoluci\363n del error global (diferencia entre soluci\363n" }} {PARA 0 "" 0 "" {TEXT -1 94 "verdadera y la aproximaci\363n) frente a \+ la del paso h podremos decidir cual es esa \"velocidad\" " }}{PARA 0 "" 0 "" {TEXT -1 24 "para cada uno de ellos. " }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 90 "La teor\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 1 " " }}{PARA 0 "" 0 "" {TEXT -1 90 "donde p es el orden del m\351todo. Por ejemplo, u n m\351todo de orden 1 indica que si reducimos " }}{PARA 0 "" 0 "" {XPPEDIT 18 0 "h;" "6#%\"hG" }{TEXT -1 91 " a la mitad entonces reduc imos el error a la mitad, y en un m\351todo de orden 2 con la misma" } }{PARA 0 "" 0 "" {TEXT -1 64 "reduci\363n de h obtenemos una reduci \363n de 4 veces el error global." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }} {PARA 0 "" 0 "" {TEXT -1 95 "En este ejercicio vamos a determinar gr \341ficamente el orden de varios m\351todos \"inc\363gnita\". Para " } }{PARA 0 "" 0 "" {TEXT -1 89 "ello vamos a correr dichos m\351todos co n la ecuaci\363n de siempre y'(t)=y(t), y(0)=1.0, cuya " }}{PARA 0 "" 0 "" {TEXT -1 92 "soluci\363n es y(t)=exp(t). Calcularemos el valor qu e nos d\341 el m\351todo en t=1 para una bateria" }}{PARA 0 "" 0 "" {TEXT -1 84 "de diferentes h's y lo guardaremos en una tabla. Posterio rmente compararemos dichos " }}{PARA 0 "" 0 "" {TEXT -1 41 "resultados con el valor exacto y(1)=e . " }}{PARA 0 "" 0 "" {TEXT -1 98 "Los h' s los construimos a partir de un h inicial, dividiendolo sucesivamente por potencias de 2: " }}{PARA 0 "" 0 "" {TEXT -1 1 " " }{XPPEDIT 18 0 "h;" "6#%\"hG" }{TEXT -1 4 " , " }{XPPEDIT 18 0 "h/2;" "6#*&%\"hG\" \"\"\"\"#!\"\"" }{TEXT -1 5 " , " }{XPPEDIT 18 0 "h/(2^2);" "6#*&%\" hG\"\"\"*$\"\"#\"\"#!\"\"" }{TEXT -1 50 " , etc. \+ " }}{PARA 0 "" 0 "" {TEXT -1 89 "Para intentar es timar el orden del m\351todo representaremos los errores cometidos fre nte a " }}{PARA 0 "" 0 "" {TEXT -1 23 "los h correspondientes." }} {PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 45 "Represent aci\363n gr\341fica (gr\341fica logar\355tmica):" }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 46 "Como la relaci\363n entre el error y h tiende a :" }}{PARA 0 "" 0 "" {TEXT -1 50 " \+ e_n = K " }{XPPEDIT 18 0 "h^p;" "6#)%\"hG %\"pG" }{TEXT -1 3 " " }}{PARA 0 "" 0 "" {TEXT -1 3 " " }}{PARA 0 "" 0 "" {TEXT -1 90 "lo mejor es representar dicha gr\341fica en ejes \+ logar\355tmicos. En efecto, tomando logar\355tmos:" }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 63 " log( e_n ) = log(K) + p log ( h ) => " }}{PARA 0 "" 0 "" {TEXT -1 101 " y = a + b x, donde: \+ y=log(e_n), a=log(K), x=log(h)" }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 101 "Luego, si en el eje de las x\264s repr esentamos log(h) y en el eje de las y\264s representamos log(e_n) la \+ " }}{PARA 0 "" 0 "" {TEXT -1 87 "gr\341fica resultante ser\341 una rec ta, cuya pendiente b coincide con el orden p del m\351todo." }}{PARA 0 "" 0 "" {TEXT -1 2 " " }}{PARA 0 "" 0 "" {TEXT -1 79 "En primer lug ar definimos la ecuaci\363n diferencial y creamos 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 condicio nes del problema" }}{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 "# Creamos una sucesi\363n de N valores de h, cada uno F veces menor" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 20 "F:=2: N:=6: 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=ar ray(1..N):" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 18 "y_met=array(1..N):" } }{PARA 6 "" 1 "" {TEXT -1 0 "" }}{PARA 6 "" 1 "" {TEXT -1 50 " \+ H maximo 1.00000 ---> H m\355nimo .03125" }}}{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 0 "" }}{PARA 0 "" 0 " " {TEXT -1 89 "Ahora definimos que m\351todo vamos a usar y pintamos e l error cometido para diferentes h's." }}{PARA 0 "" 0 "" {TEXT -1 0 " " }}{PARA 0 "" 0 "" {TEXT -1 289 "Nota: Para medir la pendiente sobr e la gr\341fica (orden del m\351todo) basta pinchar con el rat\363n en un par de puntos de la zona recta, las coordenadas se muestran en la \+ parte superior izquierda de la pantalla, anotar el valor de sus cooden adas y calcular una estimaci\363n de la pendiente como " }{XPPEDIT 18 0 "Delta;" "6#%&DeltaG" }{TEXT -1 4 "y / " }{XPPEDIT 18 0 "Delta;" "6# %&DeltaG" }{TEXT -1 3 "x ." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 111 " (Indicaci\363n: para facilitar las cuentas pro bad a seleccionar dos puntos separados por 1.0 en el eje de las X " }} {PARA 0 "" 0 "" {TEXT -1 84 "y entonces la pendiente a estimar ser\341 simplemente la variaci\363n en el eje de las Y)." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 36 "# Definimos \+ que m\351todo vamos a usar " }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 12 "meto do:=rk3:" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 18 "# Precisi\363n simple" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 26 "Digits:=18: printlevel:=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],me todo):" }}{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,cold im(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:=evalf(exp(1.0),24): " }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 3 " \+ " }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 59 "# Discrepancia entre los valore s 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 81 "loglogplot(data,axes=boxed,color=red,labels=[`Interva los: h_n`,`Errores: e_n` ]);" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }} {PARA 0 "" 0 "" {TEXT -1 1 " " }{MPLTEXT 1 0 0 "" }}{PARA 13 "" 0 "" {TEXT -1 0 "" }}{PARA 13 "" 1 "" {GLPLOT2D 400 300 300 {PLOTDATA 2 "6' -%'CURVESG6#7(7$\"\"!$!3\"H/')pqAsG\"!#<7$$!3&>\")Rm&**H5I!#=$!3Hi0(*> bh@?F+7$$!3\"RizK\"**f?gF/$!3feU:(Rv'QGF+7$$!3(eV>*p)**3.*F/$!3\"*4*=l .R&)p$F+7$$!3yCfl#)*>T?\"F+$!3T`!G9qq*zXF+7$$!3)f!*>$y*\\^]\"F+$!3a$\\ >Fl=AZ&F+-%*AXESSTYLEG6#%$BOXG-%+AXESLABELSG6$%0Intervalos:~h_nG%-Erro res:~e_nG-%'COLOURG6&%$RGBG$\"*++++\"!\")F(F(-%*AXESTICKSG6$7=/$!3++++ ++++?F+Q&.1e-16\"/$!3\")=gL/+(*)p\"F+Q!6\"/$!3dP.GX(yG_\"F+F\\o/$!3hP? n3+%zR\"F+F\\o/$!3?\")Rm&**H5I\"F+Q&.5e-1Fhn/$!3Pcjh\\([=A\"F+F\\o/$!3 !\\:\"F+F\\o/$!3Uc!3I,5p4\"F+F\\o/$!38v1c!\\dd/\"F+F\\o/!\"\"Q# .1Fhn/$!3/)=gL/+(*)pF/F\\o/$!3ivL!GX(yG_F/F\\o/$!35w.s'3+%zRF/F\\o/$!3 '>\")Rm&**H5IF/Q#.5Fhn/$!3ojN;'\\([=AF/F\\o/$!3qJu&)*f>!\\:F/F\\o/$!3Y Tc!3I,5p*!#>F\\o/$!3a7v1c!\\dd%F]rF\\o/F(Q#1.Fhn/$\"3'>\")Rm&**H5IF/F \\o/$\"3QCm>ZD@rZF/F\\o/$\"3#RizK\"**f?gF/F\\o/$\"3/)=gL/+(*)pF/Q#5.Fh n/$\"3MOk$Q]7:y(F/F\\o/$\"3MoD9+/)4X)F/F\\o/$\"3)eV>*p)**3.*F/F\\o/$\" 3v[KR%4DCa*F/F\\o7X/$!3,+++++++gF+Q&1e-06Fhn/$!3#)=gL/+(*)p&F+F\\o/$!3 dP.GX(yG_&F+F\\o/$!3jP?n3+%zR&F+F\\o/$!3>\")Rm&**H5I&F+F\\o/$!3Qcjh\\( [=A&F+F\\o/$!3=Vd)*f>!\\:&F+F\\o/$!3Vc!3I,5p4&F+F\\o/$!37v1c!\\dd/&F+F \\o/$!3++++++++]F+Q&1e-05Fhn/$!3\")=gL/+(*)p%F+F\\o/$!3cP.GX(yG_%F+F\\ o/$!3iP?n3+%zR%F+F\\o/$!3?\")Rm&**H5I%F+F\\o/$!3Qcjh\\([=A%F+F\\o/$!3= Vd)*f>!\\:%F+F\\o/$!3Uc!3I,5p4%F+F\\o/$!38v1c!\\dd/%F+F\\o/$!3,+++++++ SF+Q&.1e-3Fhn/$!3\")=gL/+(*)p$F+F\\o/$!3dP.GX(yG_$F+F\\o/$!3iP?n3+%zR$ F+F\\o/$!3?\")Rm&**H5I$F+F\\o/$!3Qcjh\\([=A$F+F\\o/$!3!\\:$F+F \\o/$!3Uc!3I,5p4$F+F\\o/$!38v1c!\\dd/$F+F\\o/$!3++++++++IF+Q&.1e-2Fhn/ $!3\")=gL/+(*)p#F+F\\o/$!3dP.GX(yG_#F+F\\o/$!3hP?n3+%zR#F+F\\o/$!3?\") Rm&**H5I#F+F\\o/$!3Pcjh\\([=A#F+F\\o/$!3!\\:#F+F\\o/$!3Uc!3I,5p 4#F+F\\o/$!38v1c!\\dd/#F+F\\oFZFinF^oFao/FeoF\\oFhoF[pF^pFapFdpFgpFjpF ]q/FaqF\\oFdqFgqFjqF^r" 1 2 0 1 0 2 6 1 2 2 1.000000 45.000000 45.000000 0 }}}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}{PARA 0 "" 0 "" {TEXT -1 0 "" }{TEXT 256 9 "Ejercicio" }{TEXT -1 2 ": " }} {PARA 0 "" 0 "" {TEXT -1 435 "Rellenar los valores de la siguiente tab la con el fin de obtener el orden (p) de los tres m\351todos propuesto s (Euler, Heun y RK3). Para ello basta pinchar con el rat\363n en un p ar de puntos de las correspondientes gr\341ficas y Maple mostrar\341 l as coordenadas de esos puntos en la parte superior izquierda de la pan talla (x1,y1,x2,y2). A continuaci\363n calcular los incrementos y una \+ estimaci\363n del orden. \277De qu\351 orden resulta ser cada m\351tod o?" }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 71 " \+ x1 y1 x2 y2 " }{XPPEDIT 18 0 "Delta;" "6#%&DeltaG" }{TEXT -1 7 "x " }{XPPEDIT 18 0 "Delta;" "6#%&DeltaG" }{TEXT -1 11 "y m= " }{XPPEDIT 18 0 " Delta;" "6#%&DeltaG" }{TEXT -1 4 "y / " }{XPPEDIT 18 0 "Delta;" "6#%&D eltaG" }{TEXT -1 14 "x p " }}{PARA 0 "" 0 "" {TEXT -1 94 "-- ---------------------------------------------------------------------- ----------------------" }}{PARA 0 "" 0 "" {TEXT -1 23 "foreuler -> \+ " }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 20 "heunform -> " }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 " " 0 "" {TEXT -1 29 "rk3 -> " }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 2 ". " }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 258 "" 0 "" {TEXT -1 0 "" }{TEXT 256 30 "Compara ci\363n de varios m\351todos." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }} {PARA 0 "" 0 "" {TEXT -1 90 "En este apartado vamos a juntar en una \+ \372nica gr\341fica el comportamiento de varios m\351todos. " }}{PARA 0 "" 0 "" {TEXT -1 53 "Partimos de la misma ecuaci\363n diferencial de siempre:" }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 50 "ecu:= 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:" }}} {EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}{PARA 0 "" 0 "" {TEXT -1 0 "" }{MPLTEXT 1 0 0 "" }{TEXT -1 0 "" }{MPLTEXT 1 0 0 "" }}}{PARA 0 " " 0 "" {TEXT -1 94 "A continuaci\363n generamos un gr\341fico como el \+ anterior, pero para varios m\351todos a la vez, y con" }}{PARA 0 "" 0 "" {TEXT -1 28 "una precisi\363n 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 "# Enum eramos los m\351todos a comparar " }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 37 "met:=\{foreuler,rk2,rk3,rk4,ad ambash\}:" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 36 "col:=\{red,green,blue, orange,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 6 "N:=6: " }}{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 mucha pr ecisi\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\351tod os" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 28 "for m from 1 to N_metodos do " }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 44 "# printf(\"**** Metodo %s *****\\n\",met[m]); " }}{PARA 0 "> \+ " 0 "" {MPLTEXT 1 0 54 "# printf(\" k Paso(h) e[h] Primer Digito mal\\n\");" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 29 "# barremos las difer entes 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,coldim(a)]),24);" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 84 "# 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 " da ta:=[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 45 "displa y(seq(g[k],k=1..N_metodos),axes=boxed);" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#<'%$rk3G%)foreulerG%)adambashG%$rk4G%$rk2G" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#<'%&greenG%'orangeG%$redG%%blueG%'yellowG" }}{PARA 13 " " 1 "" {GLPLOT2D 400 300 300 {PLOTDATA 2 "6*-%'CURVESG6$7(7$\"\"!$!3&H /')pqAsG\"!#<7$$!3&>\")Rm&**H5I!#=$!3]i0(*>bh@?F+7$$!3!RizK\"**f?gF/$! 3)*fU:(Rv'QGF+7$$!3'eV>*p)**3.*F/$!3(*>*=l.R&)p$F+7$$!3yCfl#)*>T?\"F+$ !3**H\"G9qq*zXF+7$$!3)f!*>$y*\\^]\"F+$!3D\"4?Fl=AZ&F+-%'COLOURG6&%$RGB GF($\"*++++\"!\")F(-F$6$7(7$F($!3:O&*z/70P9F/7$F-$!3FduiXp#\\H$F/7$F3$ !3]*ps3\\`rd&F/7$F8$!3)Q'R9U!yt;)F/7$F=$!3?9fg9h*\\4\"F+7$FB$!3?)z-jCP TQ\"F+-FG6&FI$\")+++!)FL$\")Vyg>FLF`o-F$6$7(7$F($!3_YlGVF`tOF+7$F-$!3# z=H?m!p<[F+7$F3$!3Z_y3.6u\"*fF+7$F8$!33U\"p%3R!3=(F+7$F=$!3fT#eFo*Qx$) F+7$FB$!30$zvMGTxd*F+-FG6&FIFJF(F(-F$6$7(7$F($!3+vK]3EC-?F+7$F-$!3<\\I t(e#*)GIF+7$F3$!3H%H&e+gLVTF+7$F8$!3S,X[z#=CI&F+7$F=$!3IUZ0N$pR['F+7$F B$!3&)=[J;AzwwF+-FG6&FIF(F(FJ-F$6$7(7$F($!3E/`jpT#)4mF/7$F-$!3@7i@].#) 46F+7$F3$!3)p]6UH*HI;F+7$F8$!3+p6)GNu5>#F+7$F=$!3[X1:S)[Dx#F+7$FB$!3Ai ]n\")4OkLF+-FG6&FIFJFJF(-%*AXESSTYLEG6#%$BOXG-%+AXESLABELSG6$%!GF^t-%% VIEWG6$%(DEFAULTGFbt" 1 2 0 1 0 2 9 1 2 2 1.000000 45.000000 45.000000 0 }}}}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }{TEXT -1 0 "" }}}{PARA 0 "" 0 "" {TEXT -1 0 "" } }{PARA 0 "" 0 "" {TEXT -1 0 "" }{MPLTEXT 1 0 0 "" }{TEXT 257 10 "Ejerc icio:" }{TEXT 258 0 "" }{TEXT -1 1 " " }}{PARA 0 "" 0 "" {TEXT -1 82 " a).- Calcular el orden de cada uno de los m\351todos, utilizando la gr \341fica anterior." }}{PARA 0 "" 0 "" {TEXT -1 143 "b).- Probar a camb iar la precisi\363n de los c\341lculos (cambiando Digits, p.e. 14, 10, 6) y estudiar cual es el l\355mite inferior \"fiable\" del error." }} {PARA 0 "" 0 "" {TEXT -1 159 "c).- Probar a cambiar el rango de variac i\363n de h (cambiando N, p.e. 8, 10) y estudiar cual es el rango de h donde la gr\341fica de la funci\363n se mantiene \"recta\"." }}} {PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{SECT 0 {PARA 4 "" 0 "" {TEXT -1 1 "3" } {TEXT 256 13 ".- Ejercicio." }}}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 66 "Resolver el problema computacional n\272 17 de \+ la Hoja de problemas. " }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 56 "Para ello se puede utilizar la siguiente implementac i\363n:" }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT 269 39 "Esquema de Runge-Kutta de cuarto orden." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 80 "Definimos la funci\363n que det ermina la ecuaci\363n diferencial: (p.e. f(x,y)=2xy )." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 51 "funcion :=proc(x,y) RETURN(2*x*y) end:funcion(x,y);" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#,$*&%\"xG\"\"\"%\"yGF&\"\"#" }}}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 114 "Definimos los l\355mites del i ntervalo [x0,x1], el tama\361o de paso h , la condici\363n inicial y0 y calculamos la malla:" }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 49 "x0:=0;x1:=2;h:=1/20;N:=(x1-x0)/h;y0 :=1;y1[1]:=y0;" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#>%#x0G\"\"!" }} {PARA 11 "" 1 "" {XPPMATH 20 "6#>%#x1G\"\"#" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#>%\"hG#\"\"\"\"#?" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#>% \"NG\"#S" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#>%#y0G\"\"\"" }}{PARA 11 " " 1 "" {XPPMATH 20 "6#>&%#y1G6#\"\"\"F'" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 39 "for n from 1 to N do x1[n]:=x0+n*h; od:" }}}{PARA 0 " " 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 61 "Definimos las 4 fu nciones de evaluaci\363n: k1(x,y),...,k4(x,y):" }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 39 "k1:=proc(x,y ) RETURN(funcion(x,y)) end:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 55 "k2:=proc(x,y) RETURN(funcion(x+h/2,y+h*k1(x,y)/2)) end:" }}} {EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 55 "k3:=proc(x,y) RETURN(funcion (x+h/2,y+h*k2(x,y)/2)) end:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 51 "k4:=proc(x,y) RETURN(funcion(x+h,y+h*k3(x,y))) end:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 32 "k1(x,y);k2(x,y);k3(x,y);k4(x,y);" } }{PARA 11 "" 1 "" {XPPMATH 20 "6#,$*&%\"xG\"\"\"%\"yGF&\"\"#" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#,$*&,&%\"xG\"\"\"#F'\"#SF'F',&%\"yGF'*&F&F'F +F'#F'\"#?F'\"\"#" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#,$*&,&%\"xG\"\"\" #F'\"#SF'\"\"\",&%\"yGF'*&F%F',&F,F'*&F&F'F,F'#F'\"#?F'F0F'\"\"#" }} {PARA 11 "" 1 "" {XPPMATH 20 "6#,$*&,&%\"xG\"\"\"#F'\"#?F'F',&%\"yGF'* &,&F&F'#F'\"#SF'\"\"\",&F+F'*&F-F',&F+F'*&F&F'F+F'F(F'F(F'#F'\"#5F'\" \"#" }}}{PARA 11 "" 1 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 20 " Definimos el m\351todo:" }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 125 "for n from 1 to N-1 do y1[n+1]:=ev alf(y1[n]+(h/6)*(k1(x1[n],y1[n])+2*k2(x1[n],y1[n])+2*k3(x1[n],y1[n])+k 4(x1[n],y1[n]))); od:" }}}{PARA 0 "" 0 "" {TEXT -1 73 "Comparamos los \+ resultados: y1[40] es la soluci\363n aproximada en x1=2 y " } {XPPEDIT 18 0 "exp(x^2);" "6#-%$expG6#*$%\"xG\"\"#" }{TEXT -1 65 " es \+ la soluci\363n exacta del problema correspondiente a f(x,y)=2xy." }} {PARA 0 "" 0 "" {TEXT -1 0 "" }}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 26 "evalf(exp(2^2),15);y1[40];" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#$\"0 U9L+:)fa!#8" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#$\"+jz4Ya!\")" }}} {EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 91 "for n from 1 to N do vect[n] :=[x1[n],y1[n]]; od: plot([seq(vect[n],n=1..N)],style=[point]);" }} {PARA 13 "" 1 "" {GLPLOT2D 400 300 300 {PLOTDATA 2 "6%-%'CURVESG6%7J7$ $\"1+++++++]!#<$\"\"\"\"\"!7$$\"1+++++++5!#;$\"1+++&>Gv+\"!#:7$$\"1+++ ++++:F1$\"1+++S8??5F47$$\"1+++++++?F1$\"1+++(*>@Q5F47$$\"1+++++++DF1$ \"1+++Yl$=1\"F47$$\"1+++++++IF1$\"1+++kAW\"4\"F47$$\"1+++++++NF1$\"1++ +^o\\F6F47$$\"1+++++++SF1$\"1+++c2eq6F47$$\"1+++++++XF1$\"1+++aFS@7F47 $$F)F1$\"1+++b$>3G\"F47$$\"1+++++++bF1$\"1+++&ze)\\8F47$$\"1+++++++gF1 $\"1+++X0vH9F47$$\"1+++++++lF1$\"1+++?:'>_\"F47$$\"1+++++++qF1$\"1+++o /CG;F47$$\"1+++++++vF1$\"1+++3Cn]Mu$F47$$ \"1+++++++7F4$\"1+++[R:5UF47$$\"1++++++]7F4$\"1+++-l\")eZF47$$\"1+++++ ++8F4$\"1+++_@%fS&F47$$\"1++++++]8F4$\"1+++6)[=<'F47$$\"1+++++++9F4$\" 1+++Qse\"3(F47$$\"1++++++]9F4$\"1+++`0:m\")F47$$F7F4$\"1,++o(=SY*F47$$ \"1++++++]:F4$\"1+++zPJ-6!#97$$\"1+++++++;F4$\"1+++yjM!H\"F\\u7$$\"1++ ++++];F4$\"1+++mX-=:F\\u7$$\"1+++++++F4$\"1+++#)GL(o$F\\u7$$\"1++++++]>F4$\" 1+++)ee+Z%F\\u7$$\"\"#F-$\"1+++jz4YaF\\u-%'COLOURG6&%$RGBG$\"#5!\"\"F- F--%&STYLEG6#%&POINTG-%+AXESLABELSG6$%!GFhx-%%VIEWG6$%(DEFAULTGF\\y" 1 2 0 1 0 2 9 1 4 2 1.000000 45.000000 45.000000 0 }}}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}{PARA 11 "" 1 "" {TEXT -1 0 "" }}} {MARK "2 0" 0 }{VIEWOPTS 1 1 0 1 1 1803 }