CV4 - Gausova eliminační  metoda

Program pro výpočet řešení soustavy lineárních rovnic Gaussovou metodou

Gaussova metoda bez výběru hlavního prvku

Gauss1[ a_List,b_List]:=Module[{ i,j,k,m,ir,is,irr,amax,p,x,ap,bp},
{ m=1;EPSCH=0.00001;n=Length[a];
x=Table[0,{i,1,n}];
ap=a;
bp=b;
For[i=1,i<n,++i,    
	  {
				For[j=i,j<n	,++j,
						{
						 p=ap[[j+1,i]]/ap[[i,i]];
	           ap[[j+1]]=ap[[j+1]]-p*ap[[i]];
	           ap[[j+1,i]]=0.0;
             bp[[j+1]]=bp[[j+1]]-p*bp[[i]];
						}];
    }];
If[Abs[ap[[n,n]]]<EPSCH,{Print[ap];Print["Matice je singularni"]},
	    { For[i=n,i>1,--i,    
	           {
							 For[j=i,j>1,--j,
								{
							  p=ap[[j-1,i]]/ap[[i,i]];
	               ap[[j-1]]=ap[[j-1]]-p*ap[[i]];
	               ap[[j-1,i]]=0.0;
                 bp[[j-1]]=bp[[j-1]]-p*bp[[i]];
								}];
	           }];
   For[i=1,i<=n,++i,
              {
               x[[i]]=bp[[i]]/(ap[[i,i]]);
               }]; 
	   }];
x}
]

Gaussova metoda s výběrem hlavního prvku

Gauss[ a_List,b_List]:=Module[{ i,j,k,m,ir,is,irr,amax,p,x,ap,bp},
{ m=1;EPSCH=0.00001;n=Length[a];
irr=Table[-1,{i,1,n}];
x=Table[0,{i,1,n}];
ap=a;
bp=b;
 For[k=1,k<=n,k++, ir=1; is=1;amax=0.0;
    { For[i=1,i<=n,++i,
			 	{If[irr[[i]]<0,
				     { For[j=1,j<=n,++j,
							  { p=Abs[ap[[i,j]]];
            	       If[p>amax,
	                  {
	                  ir=i;
	                  is=j;
	                 amax=p;
	                 }];
         	     }];
             }];
         }];
     If[amax>EPSCH,
       { irr[[ir]]=is;
       For[i=1,i<=n,++i,
          {  If[ (i!=ir) && (Abs[ap[[i,is]]]>EPSCH ) ,
	            { p=ap[[i,is]]/ap[[ir,is]];
	               ap[[i]]=ap[[i]]-p*ap[[ir]];
	               ap[[i,is]]=0.0;
             	  bp[[i]]=bp[[i]]-p*bp[[ir]];
	          }];
         }];
								
      },
     {
     m=0.0;
   Print["Matice je singularni"];
     }];
	}];
 For[i=1,i<=n,++i,
   {
   ir=irr[[i]];
   x[[ir]]=bp[[i]]/(ap[[i,ir]]);
   }];
x}
]

Použítí programu gauss1 , gauss a LinearSolve pro výpočet řešení soustavy lineárních rovnic

Definice matice a pravé strany

A={{1,2},{2,1}};
b={1,1};

Gaussova metoda bez výběru hlavního prvku

xx=Gauss1[A,b]

Gaussova metoda s výběrem hlavního prvku

XX=Gauss[A,b]

Vestavěná funkce LinearSolve

LinearSolve[A,b]

Opět použití peogramů gauss1, gauss a LinearSolve pro jiný příklad

A={{1,2,1},{2,1,-1},{1,1,1}};
b={1,1,2};

xx=Gauss1[A,b]

XX=Gauss1[A,b]

LinearSolve[A,b]

Ukázka příkladu se špatně podmíněnou maticí

Definice matice a pravé strany (matice je závislá na parametru e)

B[e_]={{10,7,11,5+e},{14,10,15,7},{15,8,19,9},{15,7,19,10}}
b={1,1,1,1}

Vypočteme si číslo podmíněnosti

Norma[AA_List,n_]:=Max[Table[Sum[Abs[AA[[i,j]]],{i,1,n}],{j,1,n}]]

Norma[B[0],4]

Norma[Inverse[B[0]],4]

Kapa=Norma[B[0],4]*Norma[Inverse[B[0]],4]

Číslo podmíněnosti je 10688, řešení bude závislé na malé změně prvků matice

xx1=Gauss[B[0],b]

xx2=Gauss[B[0.025],b]

xx3=Gauss[B[0.01],b]

Jiný příklad:
Definice matice a pravé strany (pravá strana je závislá na parametru e)

b=.;

CC={{1,1},{1,1.0001}}
b[e_]={2,2.0001+e}

Vypočteme si číslo podmíněnosti

Norma[CC,2]

Norma[Inverse[CC],2]

Kapa=Norma[CC,2]*Norma[Inverse[CC],2]

Číslo podmíněnosti je 40004, řešení bude závislé na malé změně prvků pravé strany

xx1=Gauss[CC,b[0]]

xx1=Gauss[CC,b[0.0001]]

Jiný příklad:
Ukážeme si, že podmíněnost matice nazavisí na determinantu matice.
Definice matice a pravé strany (pravá strana je závislá na parametru e).

DD[e_,n_]:=Module[{},M=Table[0,{i,1,n},{j,1,n}];For[i=1,i<=n,i++,{M[[i,i]]=1;
			   For[j=i+1,j<=n,j++,M[[i,j]]=-1]}]M[[1,n]]=M[[1,n]]+e;M]
b[e_,n_]:=Table[1+e,{i,1,n}]

n=20

A[e_]=DD[e,20];
k[e_]=b[e,20];

Vypočteme si číslo podmíněnosti

Norma[A[0],20]

Norma[Inverse[A[0]],20]

Kapa=Norma[A[0],20]*Norma[Inverse[A[0]],20]

Det[A[0]]

Číslo podmíněnosti je 10485760, ale determinat matice je roven 1, řešení bude závislé na malé změně prvků matice i pravé strany

xx1=Gauss[A[0],k[0]]

xx2=Gauss[A[0.0001],k[0]]

xx2=Gauss[A[0.0001],k[0.0001]]

Jiný příklad:
Naopak někdy je determinant velký a číslo podmíněnosti malé

DD[e_,n_]:=Module[{},M=Table[0,{i,1,n},{j,1,n}];For[i=1,i<=n,i++,M[[i,i]]=10];
			 M[[1,n]]=M[[1,n]]+e;M]
b[e_,n_]:=Table[1+e,{i,1,n}]

Vypočteme si číslo podmíněnosti a determinant

A[e_]=DD[e,5];
k[e_]=b[e,5];

Norma[A[0],5]

Norma[Inverse[A[0]],5]

Kapa=Norma[A[0],5]*Norma[Inverse[A[0]],5]

Det[A[0]]

Číslo podmíněnosti je 1, ale determinat matice je roven 1, řešení nebude závislé na malé změně prvků matice i pravé strany

xx1=Gauss[A[0],k[0]]

xx2=Gauss[A[0.0001],k[0]]

xx2=Gauss[A[0.0001],k[0.0001]]