¡Descarga programacion ejercicios y más Ejercicios en PDF de Programación C solo en Docsity!
Métodos Numéricos: Resolución de ejercicios con Maple
Tema 5: Resolución aproximada de ecuaciones
Prof. Francisco Palacios EPSEM-UPC Abril 2008. Versión 1.
Ejercicio 1.
(a) Representamos conjuntamente las funciones de los dos miembros de la ecuación x=exp(-x). La raíz es la x del punto de corte de las gráficas. > plot([x,exp(-x)],x=0..1);
0.2 0.4 x 0.6 0.8 1
1
0
(c) Si exigimos 4 decimales exactos en la aproximación mediante el método de la bisección, resulta la siguiente inecuación: > inec:=1/2^n<0.510^(-4);*
inec := <
2 n^
Resolvemos graficamente la inecuación, vemos que hacen falta 15 pasos. > plot(inec,n=0..20);
(^510) n 15 20
1
0
Programa elemental para calcular 5 pasos de bisección. > f:=x->x-exp(-x); a:=0; b:=1; n:=5; for i from 1 to n do ******** iteración,i,**************; c:=evalf((a+b)/2); if evalf(f(c)f(a))<0 then b:=c else a:=c fi od;* f := x → x − e (^ − x ) a := 0 b := 1 n := 5 ******** iteración , 1 , ************** c :=.
******** iteración Page , 2 , ************** 1
c :=. ******** iteración , 3 , ************** c :=. ******** iteración , 4 , ************** c :=. ******** iteración , 5 , ************** c :=.
Ejercicio 2. Si realizamos el gráfico con Maple, podemos representar directamente y=x-exp(-x) y observar el corte con OX. Fíjate que cuando realizamos el esquema manualmente, es preferible representar por separado las curvas y=x, y=exp(-x), y observar el punto de corte. > plot(x-exp(-x),x=0..1);
0.2 0.4^ x^ 0.6 0.8 1
0 -0. -0. -0. -0.
De la gráfica, podemos estimar el valor aproximado 0.57 para la solución. Determinamos el número de iteraciones para obtener 7 decimales exactos. 'Limpiamos' la variable n que hemos utilizado en el ejercicio anterior. > n:='n';
n := n > ineq:=1/2^n<0.5E-7;
ineq := <
2 n^
> plot(ineq,n=20..30);
(^202224) n 26 28 30
1
0
Necesitamos 25 iteraciones. > n7:=25;
n7 := 25 Usamos el programa anterior para calcular el valor de las 25 iteraciones. Observa que hemos puesto : después de od , eso evita que se impriman los cálculos realizados dentro del bucle. Hemos añadido una comando print al final que imprime el valor obtenido en la última iteración. > f:=x->x-exp(-x); a:=0; b:=1; n:=n7; for i from 1 to n do ******** iteración,i,**************; c:=evalf((a+b)/2); if evalf(f(c)f(a))<0 then b:=c else a:=c fi od: print(Valor de la aproximación =,c);* f := x → x − e (^ − x )
Page 2
> f:=x->1/x-ln(x); a:=1; b:=2; n:=4; for i from 1 to n do ******** iteración,i,**************; c:=evalf((a+b)/2); if evalf(f(c)f(a))<0 then b:=c else a:=c fi od;*
f := x → −
x
ln( x )
a := 1 b := 2 n := 4 ******** iteración , 1 , ************** c :=1. ******** iteración , 2 , ************** c :=1. ******** iteración , 3 , ************** c :=1. ******** iteración , 4 , ************** c :=1. Valor de la iteración 18. > f:=x->1/x-ln(x); a:=1; b:=2; n:=18; for i from 1 to n do ******** iteración,i,**************; c:=evalf((a+b)/2); if evalf(f(c)f(a))<0 then b:=c else a:=c fi od: print(Valor de la aproximación =,c);*
f := x → −
x
ln( x )
a := 1 b := 2 n := 18 Valor de la aproximación = , 1. Resolución con fsolve. > s:=fsolve(f(x),x=a..b);
s :=1. Error, tomando como valor exacto el proporcionado por fsolve. > error:=abs(s-c);
error :=.2723 10 - Vemos que el error real es inferior a 0.5x10^(-5). También en este caso, Maple obtiene una solución exacta con solve , usando la función especial LambertW. > s:=solve(f(x),x);
s :=
LambertW 1( ) > sf:=evalf(s);
sf :=1.
Problema 4 Representación conjunta. > plot([ln(x),exp(-x)],x=1..2,color=[red,blue]);
Page 4
1 1.2 1.4 x 1.6 1.8 2
0
El intervalo tiene longitud 1, necesitamos 18 iteraciones. Los cálculos son los mismos que en el ejercicio anterior. Observa que la cota de error sólo depende de la longitud de intervalo y el número de decimales exactos exigido. Valor de las primeras 4 iteraciones. > f:=x->ln(x)-exp(-x); a:=1; b:=2; n:=4; for i from 1 to n do ******** iteración,i,**************; c:=evalf((a+b)/2); if evalf(f(c)f(a))<0 then b:=c else a:=c fi od;* f := x →ln( x )− e (^ − x ) a := 1 b := 2 n := 4 ******** iteración , 1 , ************** c :=1. ******** iteración , 2 , ************** c :=1. ******** iteración , 3 , ************** c :=1. ******** iteración , 4 , ************** c :=1. Valor de la iteración 18. > f:=x->ln(x)-exp(-x); a:=1; b:=2; n:=18; for i from 1 to n do ******** iteración,i,**************; c:=evalf((a+b)/2); if evalf(f(c)f(a))<0 then b:=c else a:=c fi od: print(Valor de la aproximación =,c);* f := x →ln( x )− e (^ − x ) a := 1 b := 2 n := 18 Valor de la aproximación = , 1. Valor calculado por fsolve y error. > s:=fsolve(f(x),x=a..b); error:=abs(s-c); s :=1. error :=.3425 10 - En este caso Maple no puede obtener una solución real con solve.
Page 5
c:=evalf((a+b)/2); if evalf(f(c)f(a))<0 then b:=c else a:=c fi od;*
f := v →.01962 − .000014 v 1.5^ −.0000115 v^2 a := 30 b := 40 n := 3 ******** iteración , 1 , ************** c :=35. ******** iteración , 2 , ************** c :=37. ******** iteración , 3 , ************** c :=38. Valor de la iteración 11. > f:=v->9.810.002-1.4E-5v^1.5-1.15E-5v^2; a:=30; b:=40; n:=11; for i from 1 to n do ******** iteración,i,**************; c:=evalf((a+b)/2); if evalf(f(c)f(a))<0 then b:=c else a:=c fi od;** f := v →.01962 − .000014 v 1.5^ −.0000115 v^2 a := 30 b := 40 n := 11 ******** iteración , 1 , ************** c :=35. ******** iteración , 2 , ************** c :=37. ******** iteración , 3 , ************** c :=38. ******** iteración , 4 , ************** c :=38. ******** iteración , 5 , ************** c :=37. ******** iteración , 6 , ************** c :=37. ******** iteración , 7 , ************** c :=37. ******** iteración , 8 , ************** c :=37. ******** iteración , 9 , ************** c :=37. ******** iteración , 10 , ************** c :=37. ******** iteración , 11 , ************** c :=37. Cálculo de la raíz con fsolve y verificación del error. > s:=fsolve(f(x),x=a..b); error:=abs(s-c); s :=37.
error :=Page.00467988^7
Intentamos resolver con solve. > solve(f(x)=0,x);
En este caso, solve también nos proporciona la solución.
Problema 6 Expresamos la ecuación en forma nomal f(x)=0 y definimos f(x). > f:=x->tan(0.1x)-9.2exp(-x);**
f := x →tan .1( x )−9.2 e (^ − x ) Existencia y unicidad. Observa que tan(0.1x) es continua en [3,4]. > f(3.);
-. > f(4.);
. Hay cambio de signo en [3,4]. > diff(f(x),x);
.1 + .1 tan .1( x )^2 +9.2 e (^ − x ) La derivada es positiva, por lo tanto la raíz es única. Intervalo de longitud 1, para 5 decimales necesitamos 18 iteraciones. > f:=x->tan(0.1x)-9.2exp(-x); a:=3; b:=4; n:=18; for i from 1 to n do ******** iteración,i,**************; c:=evalf((a+b)/2); if evalf(f(c)f(a))<0 then b:=c else a:=c fi od: print(Valor de la aproximación =,c);*
f := x →tan .1( x )−9.2 e (^ − x ) a := 3 b := 4 n := 18 Valor de la aproximación = , 3. Calculamos el valor de f(x), tiene que ser próximo a cero. > f(c);
.9867 10 - Comparación con el resultado de **fsolve.
s:=fsolve(f(x)=0,x=a..b); error:=abs(s-c);** s :=3. error :=.2176 10 - Intentamos resolver con solve. > solve(f(x)=0,x);
2.411843808 −11.09724577 I Vemos que solve no obtiene la solución real buscada.
Problema 7 Si resolvemos el problema de forma manual, dibujaremos las curvas correspondientes a los miembros de la ecuación. > plot([x,exp(-x)],x=0..1,color=[red,blue]);
Page 8
a:=x4-delta; b:=x4+delta; f(a); f(b);
t := 8 δ :=.5000000000 10 - a :=. b :=. -.79 10 - .78 10 - Siguiendo el razonamiento explicado en el resumen de teoría, podemos asegurar que la solución es correcta.
Problema 8 Expresamos la ecuación en forma normal y calculamos las primeras 4 iteraciones. > f:=x->tan(0.1x)-9.2exp(-x); f1:=D(f); x0:=3.5; n:=4; for i from 0 to n-1 do x.(i+1):=evalf(x.i-f(x.i)/f1(x.i)); od;**
f := x →tan .1( x )−9.2 e (^ − x ) f1 := x →.1 + .1 tan .1( x )^2 +9.2 e (^ − x ) x0 :=3. n := 4 x1 :=3. x2 :=3. x3 :=3. x4 :=3. Vemos que se repiten los 8 primeros decimales. Verificamos que la solución tiene 5 decimales exactos > t:=5; delta:=0.510^(-t); a:=x4-delta; b:=x4+delta; f(a); f(b);*
t := 5 δ :=.5000000000 10 - a :=3. b :=3. -.22670 10 - .22671 10 - Observamos que hay cambio de signo, por lo tanto x4 aproxima la solución con 5 decimales exactos. Comparamos con el valor de fsolve para verificar que esto es cierto. > s:=fsolve(f(x)=0); er:=abs(s-x4); s :=3. er := 0 Obviamente, el error 0 obtenido no quiere decir que nuestra solución es exacta. El motivo de este valor es que nuestro programa proporciona el mismo valor que fsolve.
Problema 9 4 primeras iteraciones.
> f:=v->9.810.002-1.4E-5v^1.5-1.15E-5*v^2; Page 10
f1:=D(f); x0:=30.; n:=4; for i from 0 to n-1 do x.(i+1):=evalf(x.i-f(x.i)/f1(x.i)); od;
f := v →.01962 − .000014 v 1.5^ −.0000115 v^2 f1 := v →− .0000210 v .5^ −.0000230 v x0 :=30. n := 4 x1 :=38. x2 :=37. x3 :=37. x4 :=37. Se piden 3 decimales exactos. Detenemos las iteraciones en x4 al haberse "estabilizado" los 6 primeros decimales. Verificamos la validez del resultado. > t:=3; delta:=0.510^(-t); a:=x4-delta; b:=x4+delta; f(a); f(b);*
t := 3 δ :=. a :=37. b :=37. .49844 10 - -.49845 10 - Hay cambio de signo, por lo tanto, podemos garantizar que x4 aproxima la solución con 3 decimales exactos.
Problema 10 Tomamos la ecuación x^2-41=0 y aplicamos el método de Newton-Raphson con valor incial x0=6. > f:=x->x^2-41; f1:=D(f); x0:=6.5; n:=4; for i from 0 to n-1 do x.(i+1):=evalf(x.i-f(x.i)/f1(x.i)); od;
f := x → x^2 − 41 f1 := x → 2 x x0 :=6. n := 4 x1 :=6. x2 :=6. x3 :=6. x4 :=6. Observamos que el error estimado |x4-x3| es inferior a 0.5·10^(-8). Verificamos la validez del resultado > t:=6; delta:=0.510^(-t); a:=x4-delta; b:=x4+delta; f(a); f(b);*
Page 11
f := x → −
x
c
f1 :=−
x^2
x +
x
c x^2
2 x − x^2 c Vemos que efectivamente, después de simplificar, la fórmula de recurrencia del método de Newton-Raphson no contien divisiones. Aproximación de 1/9. Tomamos valor inicial 0.1. > c:=9; x0:=0.1; n:=4; for i from 0 to n-1 do x.(i+1):=evalf(2x.i-(x.i)^2c); od;**
c := 9 x0 :=. n := 4 x1 :=. x2 :=. x3 :=. x4 :=. Error con 10 dígitos > er:=1/c-x4;
er := 0 Para 1/45, tomamos valor incial 0.02. > c:=45; x0:=0.02; n:=5; for i from 0 to n-1 do x.(i+1):=evalf(2x.i-(x.i)^2c); od;**
c := 45 x0 :=. n := 5 x1 :=. x2 :=. x3 :=. x4 :=. x5 :=. Error con 10 dígitos. > er:=1/c-x5;
er := 0 Para 1/678, tomamos valor incial 1/500=0. > c:=678; x0:=0.002; n:=5; for i from 0 to n-1 do x.(i+1):=evalf(2x.i-(x.i)^2c); od;**
c := 678 x0 :=. n := 5
x1 Page:=.001288 13
x2 :=. x3 :=. x4 :=. x5 :=. Error con 10 dígitos. > er:=1/c-x5;
er :=.1 10-
Problema 13 Definimos la función de iteración. > g:=x->(x+cos(x))/2;
g := x → +
x
cos( x )
Calculamos la derivada y la representamos. > g1:=diff(g(x),x);
g1 := −
sin( x )
> plot(g1,x=0..1);
0 0.2 0.4 x 0.6 0.8 1
Vemos que g' es positiva, por lo tanto, g(x) creciente, cuando x varia en [0,1], entonces g(x) varía en [g(0),g(1)]. > g(0);
1 2 > g(1.);
. Obtenemos g([0,1])=[0.5,0.7715], por lo tanto se cumple la "propiedad inclusiva". Vamos a calcular una cota superior de |g'(x)|. Como la derivada es positiva, la función objetivo es g'. > h:=g1; h1:=diff(h,x);
h := −
sin( x )
h1 :=−
cos( x )
La derivada de h es negativa, por lo tanto, h es decreciente y tomará el máximo en el inicio del intervalo. **>
M1:=subs(x=0,h);#podemos comprobar el valor en el gráfico M1:=evalf(M1);**
M1 := −
sin 0( )
M1 :=. Determinamos el número de iteraciones para obtener 5 decimales exactos, para ello exigimos (M1)^(n)(b-a)<0.5*10^(-5). Observa que en este caso es b-a=1 y M1=0.5. > n:='n';
n := n
Page 14
> ineq:=0.2^n<0.510^(-5); plot(ineq,n=0..15,thickness=3);* ineq :=.2 n^ <.5000000000 10 -
(^246) n 8 10 12 14
1
0
> Si tomamos como valor incial x0=0.73, vemos que el método converge en 8 iteraciones. A continuación calculamos las 8 primeras iteraciones a partir de x0=0.73 y los errores estimados en cada paso. > g:=x->(x+cos(x))/2; x0:=0.73; n:=8; for i from 0 to n-1 do ***********************; x.(i+1):=evalf(g(x.i)); error.(i+1):=abs(x.(i+1)-x.i); od;
g := x → +
x
cos( x )
x0 :=. n := 8 *********************** x1 :=. error1 :=. *********************** x2 :=. error2 :=. *********************** x3 :=. error3 :=. *********************** x4 :=. error4 :=. *********************** x5 :=. error5 :=.54589 10 - *********************** x6 :=. error6 :=.8909 10 - *********************** x7 :=. error7 :=.1453 10 - *********************** x8 :=.
error8 Page :=.238 10 16 -
Observamos que en la iteración 8, el error estimado es claramente inferior a 0.5*10^(-5). Problema 14 El problema es muy parecido al anterior. Definimos la función de iteración. > g:=x->(x+exp(-x))/2;
g := x → +
x
e (^ − x )
La función de iteración es claramente continua. Para estudiar si se cumplen las hipótesis del teorema de convergencia en la iteración de punto fijo, debemos verificar que:
- g(x) cae dentro de [0.5,1] para todo x en [0.5],
- El máximo de |g'(x)| sobre [0.5,1] es menor que 1. Empezamos por calcular la derivada de la función de iteración > g1:=diff(g(x),x);
g1 := −
e (^ − x )
Una forma de investigar el signo de g1 es representar la desigualdad g1>0. > plot(g1>0,x=-3..3,thickness=3);
-3 -2 -1 (^1) x 2 3
1
0
La desigualdad toma el valor 1 cuando es cierta, por lo tanto vemos que g'(x) es positiva en el intervalo [0.5,1]. La función de iteración es creciente en el intervalo. Calculamos los valores de g en los extremos del intervalo. > g(0.5);
. > g(1.);
. Cuando x varía en [0.5,1], g(x) toma valores en [g(0.5),g(1)]. Por lo tanto se cumple la condición (1). Ahora vamos a calcular M1=max|g'(x)|. Hemos visto que la derivada es positiva en el intervalo del problema, por lo tanto |g'(x)|=g'(x). La función objetivo es h(x)=g'(x). > h:=g1;#la derivada es positiva h1:=diff(h,x);
h := −
e (^ − x )
h1 :=
e (^ − x )
h1 es claramente positiva. > plot(h1,x=0.5..1,y=0..0.4);
0.5 0.6 0.7 x 0.8 0.9 1
y
0
Por lo tanto h es creciente y toma el máximo en x=1. > M1:=subs(x=1.,h);
Page 17
x7 :=. error7 :=. *********************** x8 :=. error8 :=.36773 10 - *********************** x9 :=. error9 :=.7958 10 - *********************** x10 :=. error10 :=.1722 10 - Comparamos el valor obtenido con el que se obtiene con fsolve. > s:=fsolve(x=g(x)); err10:=abs(s-x10); s :=. err10 :=.477 10 - Vemos que, efectivamente, la aproximación x10 tiene los 5 decimales exactos exigidos. Si observamos el valor de los errores estimados, vemos que un programa controlado por el error estimado habría detenido las iteraciones en x8. En ese caso el error real, sería
> err8:=abs(s-x8);
err8 :=.10157 10 -
Punto fijo x-Lf(x)
Ejercicio 15 En este ejercicio nos proporcionan el intervalo que contiene la solución. En primer lugar debemos expresar la ecuación en forma normal f(x)=0 e identificar f(x). El programa siguiente calcula las iteraciones para la forma de punto fijo g( x )= x −λ f( x ). Como valor incial se ha tomado el punto medio del intervalo. Para mayor claridad, observa que primero contruimos g(x) como expresión y, despues, la transformamos en función con unapply. > f:=x->tan(0.1x)-9.2exp(-x); a:=3; b:=4; fa:=evalf(f(a)); fb:=evalf(f(b)); f1:=(fb-fa)/(b-a); lambda:=1/f1; g:=x-lambdaf(x); g:=unapply(g,x); ******************; x0:=evalf((a+b)/2); n:=9; for i from 0 to n-1 do x.(i+1):=g(x.i); od;* f := x →tan .1( x )−9.2 e (^ − x ) a := 3 b := 4 fa :=-. fb :=. f1 :=. λ :=2. g := x (^) − 2.481425782 tan .1( x )+22.82911719 (^) e (^ − x ) g := x → x − 2.481425782 tan .1( x )+22.82911719 e (^ − x )
****************** Page 19
x0 :=3. n := 9 x1 :=3. x2 :=3. x3 :=3. x4 :=3. x5 :=3. x6 :=3. x7 :=3. x8 :=3. x9 :=3. En el enunciado se pide una precisión de 6 decimales exactos. Tomando un margen de seguridad, hemos detenido las iteraciones cuando se han "estabilizado" los primeros 8 decimales. El resultado con fsolve es: > s:=fsolve(f(x)=0,x=a..b);
s :=3. El error exacto es: > er9:=abs(s-x9);
er9 :=.1 10 - Ejercicio 16 Resolvemos de forma similar al ejercicio anterior. > f:=x->x-cos(x); a:=0; b:=1; fa:=evalf(f(a)); fb:=evalf(f(b)); f1:=(fb-fa)/(b-a); lambda:=1/f1; g:=x-lambdaf(x); g:=unapply(g,x); ******************; x0:=evalf((a+b)/2); n:=10; for i from 0 to n-1 do x.(i+1):=g(x.i); od;* f := x → x −cos( x ) a := 0 b := 1 fa :=-1. fb :=. f1 :=1. λ :=. g :=.3149266426 x +.6850733574 cos( x ) g := x →.3149266426 x +.6850733574 cos( x ) ****************** x0 :=. n := 10 x1 :=. x2 :=. x3 :=. x4 :=. x5 :=. x6 :=. x7 :=. x8 :=.
x9 :=Page.7390851375 20
