SOLUÇÃO EXATA DO SISTEMA X=(1,0,-1,0) ATRIBUINDO dcc.ufrj.br/~rincon/Disciplinas/Metodos Numericos... ·

  • View
    216

  • Download
    0

Embed Size (px)

Text of SOLUÇÃO EXATA DO SISTEMA X=(1,0,-1,0) ATRIBUINDO dcc.ufrj.br/~rincon/Disciplinas/Metodos...

MTODOS NUMRICOS PARA RESOLUO DE SISTEMAS LINEARES (1,0,-1,0)SOLUO EXATA DO SISTEMA X=(1,0,-1,0)

> with(linalg):> restart:> with(plots):> > a[11]*x_1+a[12]*x_2+a[13]*x_3+a[14]*x_4=b[1]:> a[21]*x_1+a[22]*x_2+a[23]*x_3+a[24]*x_4=b[2]:> a[31]*x_1+a[32]*x_2+a[33]*x_3+a[34]*x_4=b[3]:> a[41]*x_1+a[42]*x_2+a[43]*x_3+a[44]*x_4=b[4]:

ATRIBUINDO VALORES> a[11]:=6:;a[12]:=2:;a[13]:=1:;a[14]:=-2:;> a[21]:=2:;a[22]:=8:;a[23]:=-1:;a[24]:=-1:;> a[31]:=1:;a[32]:=-1:;a[33]:=5:;a[34]:=1:;> a[41]:=-2:;a[42]:=-1:;a[43]:=1:;a[44]:=6:;> b[1]:=5:;b[2]:=3:;b[3]:=-4:; b[4]:=-3:;> > >

MTODOS DIRETOS

1 - MTODO DE ELIMINAO DE GAUSS COM PIVOTEAMENTO> A:=matrix(4,4,[[a[11],a[12],a[13],a[14]],[a[21],a[22],a[23],a[24]]

,[a[31],a[32],a[33],a[34]],[a[41],a[42],a[43],a[44]]]):;> b:=vector(4,[b[1],b[2],b[3],b[4]]):;> A_0:=matrix(4,5,[[a[11],a[12],a[13],a[14],b[1]],[a[21],a[22],a[23]

,a[24],b[2]],[a[31],a[32],a[33],a[34],b[3]],[a[41],a[42],a[43],a[44],b[4]]]);

:= A_0

6 2 1 -2 52 8 -1 -1 31 -1 5 1 -4

-2 -1 1 6 -3>

A_0 a matriz aumentada do sistema linear> L_1:=vector(5,[a[11],a[12],a[13],a[14],b[1]]);

:= L_1 [ ], , , ,6 2 1 -2 5> L_2:=vector(5,[a[21],a[22],a[23],a[24],b[2]]);

:= L_2 [ ], , , ,2 8 -1 -1 3> L_3:=vector(5,[a[31],a[32],a[33],a[34],b[3]]);

:= L_3 [ ], , , ,1 -1 5 1 -4> L_4:=vector(5,[a[41],a[42],a[43],a[44],b[4]]);

:= L_4 [ ], , , ,-2 -1 1 6 -3Page 1

> m_21:=a[21]/a[11];

:= m_211

3> m_31:=a[31]/a[11];

:= m_311

6> m_41:=a[41]/a[11];

:= m_41-1

3> L_22:=evalm(L_2-m_21*L_1);

:= L_22 , , , ,022

3

-4

3

-1

3

4

3> > L_32:=evalm(L_3-m_31*L_1);

:= L_32 , , , ,0-4

3

29

6

4

3

-29

6> L_42:=evalm(L_4-m_41*L_1);

:= L_42 , , , ,0-1

3

4

3

16

3

-4

3> m32:=(-4/3)/(22/3);

:= m32-2

11> m42:=(-1/3)/(22/3);

:= m42-1

22> L_33:=evalm(L_32-m32*L_22);

:= L_33 , , , ,0 0101

22

14

11

-101

22> L_43:=evalm(L_42-m42*L_22);

:= L_43 , , , ,0 014

11

117

22

-14

11> m43:=(14/11)/(101/22);

:= m4328

101> L_44:=evalm(L_43-m43*L_33);

:= L_44 , , , ,0 0 01003

2020

> A_4:=matrix(4,5,[[6,2,1,-2,5],[0, 22/3, -4/3, -1/3, 4/3],[0, 0, 101/22, 14/11, -101/22],[0, 0, 0, 1003/202, 0]]);

Page 2

>

:= A_4

6 2 1 -2 5

022

3

-4

3

-1

3

4

3

0 0101

22

14

11

-101

22

0 0 01003

2020

> x_4:=0/1003/202;

:= x_4 0> x_3:=(-101/22)/(101/22);

:= x_3 -1> x_2:=4/3-(-4/3*x_3-1/3*x_4);

:= x_2 0> x_1:=(5-(2*x_2+1*x_3-2*x_4))/6;

:= x_1 1> x:=(x_1,x_2,x_3,x_4);

:= x , , ,1 0 -1 0SOLUO DO SISTEMA PELO MEG X=(1,0,-1,0)>

DECOMPOSIO LUNOTE QUE A MATRIZ U=A_4(Matriz obtida pelo MEG) L matriz triangular inferior com diagonal unitria, formada pelos multiplicadores do MEG a menos de trocas de linhas> > U:=

LUdecomp(A,L='l',U='u',U1='u1',R='r',P='p',det='d',rank='ran');

:= U

6 2 1 -2

022

3

-4

3

-1

3

0 0101

22

14

11

0 0 01003

202> evalm(l);

Page 3

1 0 0 01

31 0 0

1

6

-2

111 0

-1

3

-1

22

28

1011

> AA:=multiply(l,U);

:= AA

6 2 1 -22 8 -1 -11 -1 5 1

-2 -1 1 6RESOLVENDO OS SISTEMAS: Ly=b e DEPOIS Ux=yLy=b (resolvendo por substituio obtemos )

> > y:=linsolve(l,b);

:= y , , ,54

3

-101

220

> x:=linsolve(U,y);

:= x [ ], , ,1 0 -1 0

MTODO DE CROUT ou L^T DL Decomposio de uma matriz simtrica A= L^T DLL: Matriz triangular inferior. e D matriz diagonal> d[11]:=a[11];

:= d11 6

> c[21]:=a[21]/d[11]:; c[31]:=a[31]/a[11]:; c[41]:=a[41]/d[11]:;> aux[21]:=c[21]*d[11]:;> d[22]:=a[22]-c[21]*aux[21]:; c[32]:=(a[32]-c[31]*aux[21])/d[22]:;

c[42]:=(a[42]-(c[41]*aux[21]))/d[22]:;> aux[31]:=c[31]*d[11]:;aux[32]:=c[32]*d[22]:;> d[33]:=a[33]-(c[31]*aux[31]+c[32]*aux[32]):;c[43]:=(a[43]-(c[41]*a

ux[31]+c[42]*aux[32]))/d[33]:;> aux[41]:=c[41]*d[11]:;aux[42]:=c[42]*d[22]:;aux[43]:=c[43]*d[33]:;> d[44]:=a[44]-(c[41]*aux[41]+c[42]*aux[42]+c[43]*aux[43]):;> LL:=matrix(4,4,[[1,0,0,0],[c[21],1,0,0],[c[31],c[32],1,0],[c[41],c

[42],c[43],1]]);

Page 4

:= LL

1 0 0 01

31 0 0

1

6

-2

111 0

-1

3

-1

22

28

1011

> DD:=matrix(4,4,[[d[11],0,0,0],[0,d[22],0,0],[0,0,d[33],0],[0,0,0,d[44]]]);

:= DD

6 0 0 0

022

30 0

0 0101

220

0 0 01003

202> UU:=multiply(DD,transpose(LL));

:= UU

6 2 1 -2

022

3

-4

3

-1

3

0 0101

22

14

11

0 0 01003

202Note que a multiplicao da matriz diagonal D, vezes a transposta da matriz L d exatamente o valor a matriz triangular superior superior U,obtida da decomposio LU., ou seja UU=D*L^t=U, logo tem-se que L*UU=L*D*L^t=ARESOLVENDO O SISTEMA LINEAR: AX=b, L*D*L^t x=b, Lz=b, Dy=z e (L^t) x = y: x a soluo.NOTE QUE O det(A)=det(L*D*L^t)=det(D)=d11*d22*d33*d44, poi s det(LL)=1.A multiplicao abaixo para obter para simples conferencia que a decomposio da matriz A esta correta.> AA:=multiply(LL,UU);

:= AA

6 2 1 -22 8 -1 -11 -1 5 1

-2 -1 1 6> zz:=linsolve(LL,b);

:= zz , , ,54

3

-101

220

> yy:=linsolve(DD,zz);Page 5

:= yy , , ,5

6

2

11-1 0

> xx:=linsolve(transpose(LL),yy);

:= xx [ ], , ,1 0 -1 0MTODO DE CHOLESKY A=L*(transposta de L)Matriz dos coeficientes A simetrica e definida positiva.

> >

L:=cholesky(A);

:= L

6 0 0 01

36

1

366 0 0

1

66

2

3366

1

222222 0

1

36

1

6666

14

11112222

1

202202606

Ly=b e transposta(L)x=y> AB:=multiply(L,transpose(L));

:= AB

6 2 1 -22 8 -1 -11 -1 5 1

-2 -1 1 6> y1:=multiply(inverse(L),b);

:= y1 , , ,5

66

2

3366

1

222222 0

> xx:=multiply(inverse(transpose(L)),y1);

:= xx [ ], , ,1 0 -1 0>

MTODOS ITERATIVOS1 - MTODO DE GAUSS-JACOBI

> > 6*x_1+2*x_2+x_3-2*x_4=5:> 2*x_1+8*x_2-x_3-x_4=3:> x_1-x_2+5*x_3+x_4=-4:> -2*x_1-x_2+1*x_3+6*x_4=-3:> with(linalg):> > x1:=(5-(2*x_2+x_3-2*x_4))/6:> x2:=(3-(2*x_1-x_3-x_4))/8:

Page 6

> x3:=(-4-(x_1-x_2+x_4))/5:> x4:=(-3-(-2*x_1-x_2+1*x_3))/6:>

INICIANDO O PROCESSO ITERATIVO> x1[1]:=5:> x2[1]:=3:> x3[1]:=-4:> x4[1]:=-3:> for i from 1 to 20 do> x1[i+1]:=evalf((5-(2*x2[i]+x3[i]-2*x4[i]))/6,5):> x2[i+1]:=evalf((3-(2*x1[i]-x3[i]-x4[i]))/8,5):> x3[i+1]:=evalf((-4-(x1[i]-x2[i]+x4[i]))/5,5):> x4[i+1]:=evalf((-3-(-2*x1[i]-x2[i]+x3[i]))/6,5):> E[i]:=evalf(max(abs(x1[i+1]-x1[i]),abs(x2[i+1]-x2[i]),abs(x3[i+1]-

x3[i]),abs(x4[i+1]-x4[i])),8): > X[i+1]:=evalf(vector(4,[x1[i+1],x2[i+1],x3[i+1],x4[i+1]]),5);> od;>

:= x12 -.50000

:= x22 -1.7500

:= x32 -.60000

:= x42 2.3333

:= E1 5.50000

:= X2 [ ], , ,-.50000 -1.7500 -.60000 2.3333

:= x13 2.2945

:= x23 .71666

:= x33 -1.5167

:= x43 -.85834

:= E2 3.19164

:= X3 [ ], , ,2.2945 .71666 -1.5167 -.85834

:= x14 .56111

:= x24 -.49551

:= x34 -.94393

:= x44 .63705

:= E3 1.73339

:= X4 [ ], , ,.56111 -.49551 -.94393 .63705

Page 7

:= x15 1.3682

:= x25 .19636

:= x35 -1.1387

:= x45 -.23823

:= E4 .87528

:= X5 [ ], , ,1.3682 .19636 -1.1387 -.23823

:= x16 .87825

:= x26 -.13917

:= x36 -.98665

:= x46 .17858

:= E5 .48995

:= X6 [ ], , ,.87825 -.13917 -.98665 .17858

:= x17 1.1037

:= x27 .054433

:= x37 -1.0392

:= x47 -.06601

:= E6 .24459

:= X7 [ ], , ,1.1037 .054433 -1.0392 -.06601

:= x18 .96639

:= x28 -.039081

:= x38 -.99660

:= x48 .05017

:= E7 .13731

:= X8 [ ], , ,.96639 -.039081 -.99660 .05017

:= x19 1.0292

:= x29 .015091

:= x39 -1.0111

:= x49 -.01828

:= E8 .06845

:= X9 [ ], , ,1.0292 .015091 -1.0111 -.01828

:= x110 .99073

:= x210 -.010975Page 8

:= x310 -.99914

:= x410 .01411

:= E9 .03847

:= X10 [ ], , ,.99073 -.010975 -.99914 .01411

:= x111 1.0082

:= x211 .0041938

:= x311 -1.0031

:= x411 -.00507

:= E10 .01918

:= X11 [ ], , ,1.0082 .0041938 -1.0031 -.00507

:= x112 .99742

:= x212 -.0030738