Upload
others
View
9
Download
0
Embed Size (px)
Citation preview
Practica: Métodos iterativos para sistemas
Operaciones con Matrices y vectores
• Definición de matriz
• Definición de vector
• Matriz por vector> a:=matrix(3,3,[1,0,3,1,2,1,3,2,2]);
:= a⎡
⎣
⎢⎢⎢
⎤
⎦
⎥⎥⎥
1 0 31 2 13 2 2
> v:=[1,2,1];
:= v [ ], ,1 2 1> w=evalm(a&*v);
= w [ ], ,4 6 9
• Inversa
• Producto de matrices
• Evaluación float > b:=matrix(3,3,[1,0,1,0,1,1,1,1,0]);
:= b⎡
⎣
⎢⎢⎢
⎤
⎦
⎥⎥⎥
1 0 10 1 11 1 0
> b1:=evalm(b^(-1));
:= b1
⎡
⎣
⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢
⎤
⎦
⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥
12
-12
12
-12
12
12
12
12
-12
> c:=b&*b1;
:= c &*b b1> evalm(c);
⎡
⎣
⎢⎢⎢
⎤
⎦
⎥⎥⎥
1 0 00 1 00 0 1
> evalf(b1); #no funciona
b1> b1f:=evalf(evalm(b1));
:= b1f⎡
⎣
⎢⎢⎢
⎤
⎦
⎥⎥⎥
.5000000000 -.5000000000 .5000000000-.5000000000 .5000000000 .5000000000.5000000000 .5000000000 -.5000000000
> >
• Producto por escalar
• Suma de matrices> a:=matrix(3,3,[1,1,1,2,1,2,3,4,1]);
:= a⎡
⎣
⎢⎢⎢
⎤
⎦
⎥⎥⎥
1 1 12 1 23 4 1
> c:=3*a;
:= c 3 a> evalm(c);
Page 1
⎡
⎣
⎢⎢⎢
⎤
⎦
⎥⎥⎥
3 3 36 3 69 12 3
> b:=matrix(3,3,[1,1,2,1,1,1,0,0,1]);
:= b⎡
⎣
⎢⎢⎢
⎤
⎦
⎥⎥⎥
1 1 21 1 10 0 1
> d:=a+b;
:= d + a b> evalm(a+b);
⎡
⎣
⎢⎢⎢
⎤
⎦
⎥⎥⎥
2 2 33 2 33 4 2
• Vectores como matrices columna> a:=matrix(3,3,[1,0,1,1,1,1,2,0,0]);
:= a⎡
⎣
⎢⎢⎢
⎤
⎦
⎥⎥⎥
1 0 11 1 12 0 0
> v:=matrix(3,1,[1,2,1]);
:= v⎡
⎣
⎢⎢⎢
⎤
⎦
⎥⎥⎥
121
> w:=a&*v;
:= w &*a v> evalm(w);
⎡
⎣
⎢⎢⎢
⎤
⎦
⎥⎥⎥
242
>
Page 2
Práctica Métodos iterativos para sistemas de ecuaciones 1
Ejercicio 1Dadas las matrices
A =
⎛⎜⎝ 1 2 31 −1 12 1 −1
⎞⎟⎠ , B =
⎛⎜⎝ −1/5 1/7 01/8 −1/17 1/240 1/27 −1/121
⎞⎟⎠ ,y los vectores
c =
⎛⎜⎝ 1.22.57.1
⎞⎟⎠ , v =
⎛⎜⎝ 123
⎞⎟⎠Calcula de foma exacta y en forma de aproximación decimal(a) D = AB(b) A1 = A
−1
(c) H = 13A+ 2
3B
(d) w = Bv + c(e) Si (
u(0) = v
u(j+1) = Bu(j) + c
calcula una aproximación decimal de u(18).
Solución
(a) D =
⎛⎜⎝ .0 5 . 13632 11951 .0 58539 9449−. 325 . 23871 77093 −4. 99311 2948× 10−2−. 275 . 18985 37193 4. 99311 2948× 10−2
⎞⎟⎠(b) A1 =
⎛⎜⎝ 0 . 33333 33333 . 33333 33333. 2 −. 46666 66667 . 13333 33333. 2 . 2 −. 2
⎞⎟⎠(c) H =
⎛⎜⎝ . 2 . 76190 47619 1.0. 41666 66667 −. 37254 90196 . 36111 11111. 66666 66667 . 35802 46914 −. 33884 29752
⎞⎟⎠(d) w =
⎛⎜⎝ 1. 28571 42862. 63235 29417. 14928 0686
⎞⎟⎠(e) u(18) =
⎛⎜⎝ 1.3332939132. 7996688687. 144644788
⎞⎟⎠
Practica: Métodos iterativos para sistemas
Ejercicio 1
> a:=matrix(3,3,[1,2,3,1,-1,1,2,1,-1]);
:= a⎡
⎣
⎢⎢⎢
⎤
⎦
⎥⎥⎥
1 2 31 -1 12 1 -1
> b:=matrix(3,3,[-1/5,1/7,0,1/8,-1/17,1/24,0,1/27,-1/121]);
:= b
⎡
⎣
⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢
⎤
⎦
⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥
-15
17
0
18
-117
124
01
27-1
121> c:=[1.2,2.5,7.1];
:= c [ ], ,1.2 2.5 7.1> v:=[1,2,3];
:= v [ ], ,1 2 3> d:=evalm(a&*b);
:= d
⎡
⎣
⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢
⎤
⎦
⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥
120
1461071
851452
-1340
7673213
-1452904
-1140
6103213
1452904
> df:=evalf(evalm(d));
:= df⎡
⎣
⎢⎢⎢
⎤
⎦
⎥⎥⎥
.05000000000 .1363211951 .05853994490-.3250000000 .2387177093 -.04993112948-.2750000000 .1898537193 .04993112948
> a1:=evalm(a^(-1));
:= a1
⎡
⎣
⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢
⎤
⎦
⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥
013
13
15
-715
215
15
15
-15
> a1f:=evalf(evalm(a1));
:= a1f⎡
⎣
⎢⎢⎢
⎤
⎦
⎥⎥⎥
0 .3333333333 .3333333333.2000000000 -.4666666667 .1333333333.2000000000 .2000000000 -.2000000000
> h:=1/3*a+2/3*b;
:= h + 13
a23
b
> hf:=evalf(evalm(h));
:= hf⎡
⎣
⎢⎢⎢
⎤
⎦
⎥⎥⎥
.2000000000 .7619047619 1.
.4166666667 -.3725490196 .3611111111
.6666666667 .3580246914 -.3388429752> w:=evalm(b&*v+c);
:= w [ ], ,1.285714286 2.632352941 7.149280686> u0:=v;
n:=17; for i from 0 to n do u.(i+1):=evalf(evalm(b&*u.i+c)); od;
:= u0 [ ], ,1 2 3 := n 17Page 1
:= u1 [ ], ,1.285714286 2.632352941 7.149280686 := u2 [ ], ,1.318907563 2.803756691 7.138409589 := u3 [ ], ,1.336755158 2.797370314 7.144847720 := u4 [ ], ,1.332273299 2.800245188 7.144557980 := u5 [ ], ,1.333580367 2.799503773 7.144666851 := u6 [ ], ,1.333213037 2.799715306 7.144638491 := u7 [ ], ,1.333316722 2.799655765 7.144646560 := u8 [ ], ,1.333287479 2.799672564 7.144644289 := u9 [ ], ,1.333295728 2.799667826 7.144644929 := u10 [ ], ,1.333293401 2.799669162 7.144644749 := u11 [ ], ,1.333294057 2.799668785 7.144644800 := u12 [ ], ,1.333293872 2.799668891 7.144644785 := u13 [ ], ,1.333293924 2.799668861 7.144644789 := u14 [ ], ,1.333293910 2.799668870 7.144644788 := u15 [ ], ,1.333293914 2.799668868 7.144644788 := u16 [ ], ,1.333293913 2.799668868 7.144644788 := u17 [ ], ,1.333293913 2.799668868 7.144644788 := u18 [ ], ,1.333293913 2.799668868 7.144644788
>
Page 2
Practica: Métodos iterativos para sistemas
Libreria de álgebra lineal linalg
> with(linalg);Warning, new definition for normWarning, new definition for trace
BlockDiagonal GramSchmidt JordanBlock LUdecomp QRdecomp Wronskian addcol addrow adj adjoint angle, , , , , , , , , , ,[augment backsub band basis bezout blockmatrix charmat charpoly cholesky col coldim colspace colspan companion, , , , , , , , , , , , , ,concat cond copyinto crossprod curl definite delcols delrows det diag diverge dotprod eigenvals eigenvalues, , , , , , , , , , , , , ,eigenvectors eigenvects entermatrix equal exponential extend ffgausselim fibonacci forwardsub frobenius gausselim, , , , , , , , , , ,gaussjord geneqns genmatrix grad hadamard hermite hessian hilbert htranspose ihermite indexfunc innerprod, , , , , , , , , , , ,intbasis inverse ismith issimilar iszero jacobian jordan kernel laplacian leastsqrs linsolve matadd matrix minor, , , , , , , , , , , , , ,minpoly mulcol mulrow multiply norm normalize nullspace orthog permanent pivot potential randmatrix randvector, , , , , , , , , , , , ,rank ratform row rowdim rowspace rowspan rref scalarmul singularvals smith stack submatrix subvector sumbasis, , , , , , , , , , , , , ,swapcol swaprow sylvester toeplitz trace transpose vandermonde vecpotent vectdim vector wronskian, , , , , , , , , , ]
• determinante, inversa> a:=matrix(2,2,[1,2,4,5]);
:= a ⎡⎣⎢
⎤⎦⎥
1 24 5
> a1:=inverse(a);
:= a1
⎡
⎣
⎢⎢⎢⎢⎢⎢
⎤
⎦
⎥⎥⎥⎥⎥⎥
-53
23
43
-13
> d:=det(a);
:= d -3> evalm(d*a1);
⎡⎣⎢
⎤⎦⎥
5 -2-4 1
> b:=matrix(3,3,[1,2,x,1,x,2,0,1,0]);
:= b⎡
⎣
⎢⎢⎢
⎤
⎦
⎥⎥⎥
1 2 x1 x 20 1 0
> det(b);
− + 2 x
• resolucion de sistemas lineales> a:=matrix(3,3,[1,2,1,2,-1,2,1,3,3]);
:= a⎡
⎣
⎢⎢⎢
⎤
⎦
⎥⎥⎥
1 2 12 -1 21 3 3
> b:=[1,0,1];
:= b [ ], ,1 0 1> s:=linsolve(a,b);
:= s⎡
⎣⎢⎢
⎤
⎦⎥⎥, ,
25
25
-15
> evalm(a&*s);
[ ], ,1 0 1
• Función vectorial, jacobiano> f:=(x,y,z)->[x+y^2,x*y-z,z^2+y^2];
:= f → ( ), ,x y z [ ], , + x y2 − x y z + z2 y2
> f(1,2,3);
[ ], ,5 -1 13> jf:=jacobian(f(x,y,z),[x,y,z]);
Page 1
:= jf⎡
⎣
⎢⎢⎢
⎤
⎦
⎥⎥⎥
1 2 y 0y x -10 2 y 2 z
• Sustitución en el jacobiano> v:=[1,2,3];
:= v [ ], ,1 2 3> jfv:=subs({x=v[1],y=v[2],z=v[3]},evalm(jf));
:= jfv⎡
⎣
⎢⎢⎢
⎤
⎦
⎥⎥⎥
1 4 02 1 -10 4 6
>
Page 2
Práctica Métodos iterativos para sistemas de ecuaciones 2
Ejercicio 2: Libreria linalg
(a) Resuelve el sistema, calcula una aproximación decimal de la solución.⎧⎪⎨⎪⎩10x+ y + z = 12−x+ 25y − 2z = −232x− 3y + 24z = 16
(b) Define la función vectorial
F
⎛⎜⎝ xyz
⎞⎟⎠ =⎛⎜⎝ x3 − 2xzx− 2xy + z2x2 − yz
⎞⎟⎠
calcula la imágen del vector v =
⎛⎜⎝ 1−10
⎞⎟⎠ .(c) Calcula el jacobiano de F y su valor en v.(d) Resuelve el sistema, calcula una aproximación decimal de la solución.⎧⎪⎨⎪⎩
10x+ 7y + z = −12−5x+ 25y − 2z = −1232x− 13y + 24z = 160
Solución
(a) x = 1. 23744 2922, y = −0. 83375 61306, z = 0. 45932 69068
(b) F
⎛⎜⎝ 1−10
⎞⎟⎠ =⎛⎜⎝ 131
⎞⎟⎠(c) JF =
⎛⎜⎝ 3x2 − 2z 0 −2x1− 2y −2x 2z2x −z −y
⎞⎟⎠, JF
⎛⎜⎝ 1−10
⎞⎟⎠ =⎛⎜⎝ 3 0 −23 −2 02 0 1
⎞⎟⎠(d) x = 1. 39104 614, y = −4. 30424 8515, z = 4. 21927 8209
Practica: Métodos iterativos para sistemas
Ejercicio 2
*********** (a) *********************> with(linalg):
a:=matrix(3,3,[10,1,1,-1,25,-2,2,-3,24]); b:=[12,-23,16];
:= a⎡
⎣
⎢⎢⎢
⎤
⎦
⎥⎥⎥
10 1 1-1 25 -22 -3 24
:= b [ ], ,12 -23 16> s:=linsolve(a,b);
:= s⎡
⎣⎢⎢
⎤
⎦⎥⎥, ,
271219
-49305913
27165913
> sf:=evalf(evalm(s));
:= sf [ ], ,1.237442922 -.8337561306 .4593269068************ (b) **********************> f:=(x,y,z)->[x^3-2*x*z,x-2*x*y+z^2,x^2-y*z];
:= f → ( ), ,x y z [ ], , − x3 2 z x − + x 2 x y z2 − x2 y z> f(1,-1,0);
[ ], ,1 3 1************ (c) **********************> jf:=jacobian(f(x,y,z),[x,y,z]);
:= jf⎡
⎣
⎢⎢⎢⎢
⎤
⎦
⎥⎥⎥⎥
− 3 x2 2 z 0 −2 x − 1 2 y −2 x 2 z2 x −z −y
> jfv:=subs({x=1,y=-1,z=0},evalm(jf));
:= jfv⎡
⎣
⎢⎢⎢
⎤
⎦
⎥⎥⎥
3 0 -23 -2 02 0 1
*********** (d) *************************> a:=matrix(3,3,[10,7,1,-5,25,-2,2,-13,24]);
b:=[-12,-123,160]; s:=linsolve(a,b); sf:=evalf(evalm(s));
:= a⎡
⎣
⎢⎢⎢
⎤
⎦
⎥⎥⎥
10 7 1-5 25 -22 -13 24
:= b [ ], ,-12 -123 160
:= s⎡
⎣⎢⎢
⎤
⎦⎥⎥, ,
30452189
-94222189
92362189
:= sf [ ], ,1.391046140 -4.304248515 4.219278209>
Page 1
Práctica Métodos iterativos para sistemas de ecuaciones 3
Ejercicio 3: Método iterativo lineal
Consideramos el sistema⎧⎪⎨⎪⎩10x+ y + z = 1−x+ 25y − 2z = −132x− y + 24z = 10
(a) Define las matrices A y b de la forma matricial Ax = b(b) Define la matriz N correspondiente al método de jacobi(c) Calcula P = N−A(d) Calcula M = N−1P(e) Calcula c = N−1b(f) Estudia la convergencia(g) Determina el número de iteraciones necesarias para obtener 6 decimalesexactos, tomando x(0) = ~0(h) Calcula el valor aproximado(i) Calcula la solución exacta y verifica el resultado(j) Repite todos los pasos emplando ahora el método de Gauss-Seidel.
Solución
(a) A =
⎛⎜⎝ 10 1 1−1 25 −22 −1 24
⎞⎟⎠ , b =⎛⎜⎝ 1−1310
⎞⎟⎠ (b) N =
⎛⎜⎝ 10 0 00 25 00 0 24
⎞⎟⎠(c)P =
⎛⎜⎝ 0 −1 −11 0 2−2 1 0
⎞⎟⎠ (d)M =
⎛⎜⎝ 0 − 110− 110
125
0 225
− 112
124
0
⎞⎟⎠ (e) c =
⎛⎜⎝110
−1325512
⎞⎟⎠(f) kMk∞ = 1
5< 1, el método es convergente
(g) (0.2)j
0.8× 13
25≤ 0.5× 10−6 =⇒ j ≥ 8.747, necesitamos 9 iteraciones
(h) Calculamos con 14 decimalesx = 0.10972945715932, y = −0.48462443285167, z = 0.38732986060353(i) x = 0. 10972 945723408, y = −0. 48462 443286843, z = 0. 38732 986052764°°°α− x(9)°°° = −0.7589× 10−10(j) x(9) = [.10972945723411,−.48462443286843, .38732986052764]
Practica: Métodos iterativos para sistemas
Ejercicio 3: Método iterativo lineal
> with(linalg):Warning, new definition for normWarning, new definition for trace> a:=matrix(3,3,[10,1,1,-1,25,-2,2,-1,24]);
b:=[1,-13,10]; n:=matrix(3,3,[a[1,1],0,0,0,a[2,2],0,0,0,a[3,3]]); n1:=inverse(n); p:=evalm(n-a); m:=evalm(n1&*p); c:=evalm(n1&*b);
:= a⎡
⎣
⎢⎢⎢
⎤
⎦
⎥⎥⎥
10 1 1-1 25 -22 -1 24
:= b [ ], ,1 -13 10
:= n⎡
⎣
⎢⎢⎢
⎤
⎦
⎥⎥⎥
10 0 00 25 00 0 24
:= n1
⎡
⎣
⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢
⎤
⎦
⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥
110
0 0
0125
0
0 01
24
:= p⎡
⎣
⎢⎢⎢
⎤
⎦
⎥⎥⎥
0 -1 -11 0 2
-2 1 0
:= m
⎡
⎣
⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢
⎤
⎦
⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥
0-110
-110
125
0225
-112
124
0
:= c⎡
⎣⎢⎢
⎤
⎦⎥⎥, ,
110
-1325
512
Convergencia> nm:=norm(m,infinity);
nc:=norm(c,infinity);
:= nm15
:= nc1325
||M||<1, el método es convergente> ineq:=nm^j/(1-nm)*nc<0.5*10^(-6);
:= ineq < 1320
⎛
⎝⎜⎜
⎞
⎠⎟⎟
15
j
.5000000000 10-6
> plot(ineq,j=5..10,thickness=3);
Page 1
j1098765
1
0.8
0.6
0.4
0.2
0
gráficamente, j>8.7> Digits:=14;
a:=matrix(3,3,[10,1,1,-1,25,-2,2,-1,24]); b:=[1,-13,10]; n:=matrix(3,3,[a[1,1],0,0,0,a[2,2],0,0,0,a[3,3]]); n1:=inverse(n); p:=evalm(n-a); m:=evalm(n1&*p); c:=evalm(n1&*b); nmax:=8; x0:=[0,0,0]; for i from 0 to nmax do x.(i+1):=evalf(evalm(m&*x.i+c)); od;
>
:= Digits 14
:= a⎡
⎣
⎢⎢⎢
⎤
⎦
⎥⎥⎥
10 1 1-1 25 -22 -1 24
:= b [ ], ,1 -13 10
:= n⎡
⎣
⎢⎢⎢
⎤
⎦
⎥⎥⎥
10 0 00 25 00 0 24
:= n1
⎡
⎣
⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢
⎤
⎦
⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥
110
0 0
0125
0
0 01
24
:= p⎡
⎣
⎢⎢⎢
⎤
⎦
⎥⎥⎥
0 -1 -11 0 2
-2 1 0
:= m
⎡
⎣
⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢
⎤
⎦
⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥
0-110
-110
125
0225
-112
124
0
:= c⎡
⎣⎢⎢
⎤
⎦⎥⎥, ,
110
-1325
512
:= nmax 8Page 2
:= x0 [ ], ,0 0 0 := x1 [ ], ,.10000000000000 -.52000000000000 .41666666666667 := x2 [ ], ,.11033333333333 -.48266666666667 .38666666666667 := x3 [ ], ,.10960000000000 -.48465333333333 .38736111111112 := x4 [ ], ,.10972922222222 -.48462711111111 .38733944444445 := x5 [ ], ,.10972876666667 -.48462367555556 .38732976851852 := x6 [ ], ,.10972939070370 -.48462446785185 .38732994962963 := x7 [ ], ,.10972945182222 -.48462442840148 .38732986461420 := x8 [ ], ,.10972945637873 -.48462443275798 .38732986116476 := x9 [ ], ,.10972945715932 -.48462443285167 .38732986060353
> s:=linsolve(a,b);
:= s⎡
⎣⎢⎢
⎤
⎦⎥⎥, ,
6535951
-28845951
23055951
> sf:=evalf(evalm(s));
:= sf [ ], ,.10972945723408 -.48462443286843 .38732986052764> er9:=evalm(sf-x9);
:= er9 [ ], ,.7476 10-10 -.1676 10-10 -.7589 10-10
Método de Gauss-Seidel> Digits:=16;
a:=matrix(3,3,[10,1,1,-1,25,-2,2,-1,24]); b:=[1,-13,10]; n:=matrix(3,3,[a[1,1],0,0,a[2,1],a[2,2],0,a[3,1],a[3,2],a[3,3]]); n1:=inverse(n); p:=evalm(n-a); m:=evalm(n1&*p); c:=evalm(n1&*b); nmax:=8; x0:=[0,0,0]; for i from 0 to nmax do x.(i+1):=evalf(evalm(m&*x.i+c)); od;
:= Digits 16
:= a⎡
⎣
⎢⎢⎢
⎤
⎦
⎥⎥⎥
10 1 1-1 25 -22 -1 24
:= b [ ], ,1 -13 10
:= n⎡
⎣
⎢⎢⎢
⎤
⎦
⎥⎥⎥
10 0 0-1 25 02 -1 24
:= n1
⎡
⎣
⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢
⎤
⎦
⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥
110
0 0
1250
125
0
-496000
1600
124
:= p⎡
⎣
⎢⎢⎢
⎤
⎦
⎥⎥⎥
0 -1 -10 0 20 0 0
:= m
⎡
⎣
⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢
⎤
⎦
⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥
0-110
-110
0-1
25019250
049
600023
2000
:= c⎡
⎣⎢⎢
⎤
⎦⎥⎥, ,
110
-129250
23216000
Page 3
:= nmax 8 := x0 [ ], ,0 0 0
:= x1 [ ], ,.1000000000000000 -.5160000000000000 .3868333333333333 := x2 [ ], ,.1129166666666667 -.4845366666666667 .3870679166666666 := x3 [ ], ,.1097468750000000 -.4846446916666667 .3873275649305555 := x4 [ ], ,.1097317126736111 -.4846245262986111 .3873296686814236 := x5 [ ], ,.1097294857617188 -.4846244470750174 .3873298575583977 := x6 [ ], ,.1097294589516620 -.4846244330372617 .3873298603774756 := x7 [ ], ,.1097294572659786 -.4846244328791628 .3873298605245366 := x8 [ ], ,.1097294572354626 -.4846244328686186 .3873298605275190 := x9 [ ], ,.1097294572341100 -.4846244328684341 .3873298605276394
> er9:=evalm(sf-x9);
:= er9 [ ], ,-.300 10-13 .41 10-14 .6 10-15
>
Page 4
Práctica Métodos iterativos para sistemas de ecuaciones 4
Ejercicio 4: Método Newton-Raphson
Consideramos el sistema (x2 + y2 = 4y = (x− 1)2
(a) Representa conjutamente las curvas y calcula una aproximación de lospuntos de corte(b) Plantea el sistema como una ecación vectorial F(x) = 0(c) Define la función vectorial F(x)(d) Calcula el jacobiano(e) Escribe un programa que permita aplicar el método de Newton-Raphson,
toma como valor inicial x(0) =
Ã1.80.7
!
Practica: Métodos iterativos para sistemas
Ejercicio 4: Método de Newton-Raphson
> with(plots):> eq1:=x^2+y^2-4;
eq2:=y-(x-1)^2;
:= eq1 + − x2 y2 4
:= eq2 − y ( ) − x 1 2
> implicitplot({eq1=0,eq2=0},x=-1..4,y=-1..4);
x321-1
y
4
3
2
1
0
-1
Solución de primer cuadrante x=1.8, y=0.7;> f:=unapply([eq1,eq2],x,y);
:= f → ( ),x y [ ], + − x2 y2 4 − y ( ) − x 1 2
> with(linalg):> jf:=jacobian(f(x,y),[x,y]);
:= jf ⎡⎣⎢
⎤⎦⎥
2 x 2 y− + 2 x 2 1
> x0:=[1.8,0.7]; n:=3; for i from 0 to n do `************`,i+1,`*************`; jfx:=subs(x=x.i[1],y=x.i[2],evalm(jf)); x.(i+1):=evalm(x.i-inverse(jfx)&*f(x.i[1],x.i[2])); od;
:= x0 [ ],1.8 .7 := n 3
, ,************ 1 *************
:= jfx ⎡⎣⎢
⎤⎦⎥
3.6 1.4-1.6 1
:= x1 [ ],1.8606164383562 .73698630136986, ,************ 2 *************
:= jfx ⎡⎣⎢
⎤⎦⎥
3.7212328767124 1.4739726027397-1.7212328767124 1
:= x2 [ ],1.8589453355656 .73778429690542, ,************ 3 *************
:= jfx ⎡⎣⎢
⎤⎦⎥
3.7178906711312 1.4755685938108-1.7178906711312 1
:= x3 [ ],1.8589441280931 .73778501518413Page 1
, ,************ 4 *************
:= jfx ⎡⎣⎢
⎤⎦⎥
3.7178882561862 1.4755700303683-1.7178882561862 1
:= x4 [ ],1.8589441280924 .73778501518444>
Page 2