Sustavi linearnih jednadzbi

Potprogrami - Moduli

(* Algoritam FS *)
FS[n_,a_,b_]:=Module[{i,j},   x=Table[0,{i,n}];
		   x[[1]]=b[[1]]/a[[1,1]];
	   	Do[
		         	x[[i]]=(b[[i]]-Sum[a[[i,j]] x[[j]],{j,i-1}])/a[[i,i]]
				       ,{i,2,n}]
                                    	              ]
(* Algoritam BS *)
BS[n_,a_,b_]:=Module[{i,j,k},       x=Table[0,{i,n}];
		x[[n]]=b[[n]]/a[[n,n]];
		Do[k=n-i;
			     x[[k]]=(b[[k]]-Sum[a[[k,j]] x[[j]],{j,k+1,n}])/a[[k,k]]
    			 ,{i,1,n-1}];N[x,5]
	                                                    ]

(* LU dekompozicija - Crout's method *)
LU[a_,n_]:=Module[{i,j,k}, L=Table[0,{i,n},{j,n}];  U=Table[0,{i,n},{j,n}]; 
		Do[L[[i,i]]=1,{i,n}];
	   	Do[
                       Do[
     	U[[i,j]]=a[[i,j]]-Sum[L[[i,k]] U[[k,j]],{k,i-1}],
				{j,i,n}];
               Do[
		L[[j,i]]=(a[[j,i]]-Sum[L[[j,k]]U[[k,i]],{k,i-1}])/U[[i,i]],
					{j,i+1,n}]	
			,{i,n}];
		Print["L=",MatrixForm[L],",  U=",MatrixForm[U]]
	                                               ]

(* Cholesky dekompozicija *)
Cholesky[a_,n_]:=Module[{i,j,k,p},
L=Table[0,{i,n},{j,n}];
Do[L[[k,k]]=Sqrt[a[[k,k]]-Sum[L[[k,p]]^2,{p,k-1}]];
		Do[
		L[[i,k]]=(a[[i,k]]-Sum[L[[i,p]]L[[k,p]],{p,k-1}])/L[[k,k]],{i,k+1,n}]
	  ,{k,n}];
	Print["L=",MatrixForm[L]]
	                                                        	]

(* Rjesavanje trodijagonalnog sustava LU-dekompozicijom *)
(* b je gornja, a c donja dijagonala *)
LUtri[n_,a_,b_,c_,y_]:=Module[{d,e,f,z,x},
z=Table[0,{i,n}]; x=Table[0,{i,n}]; e=Table[0,{i,n}]; f=Table[0,{i,n}];d=Table[0,{i,n}]; 
(* Izracunavanje *)
e[[1]]=a[[1]];f[[1]]=b[[1]];d[[1]]=c[[1]]/e[[1]]; Do[e[[i]]=a[[i]]-d[[i-1]] f[[i-1]];f[[i]]=b[[i]]; d[[i]]=c[[i]]/e[[i]],{i,2,n-1}];
e[[n]]=a[[n]]-d[[n-1]] f[[n-1]];

(* Rjesavanje sustava LUx=y *)
z[[1]]=y[[1]]; 
Do[z[[i]]=y[[i]]-d[[i-1]]z[[i-1]],{i,2,n}];
x[[n]]=z[[n]]/e[[n]];
Do[x[[n-i-1]]=(z[[n-i-1]]-f[[n-i-1]] x[[n-i]] )/e[[n-i-1]],{i,0,n-2}];
x//N]

(* Loo Norma matrice a *)
norma[a_,n_]:=Module[{s},
		s=Table[0,{i,n}];
Do[s[[i]]=Sum[Abs[a[[i,j]]],{j,n}],{i,n}];
Max[s]//N
	];

(* pretvaranje broja w u k-znamenkasti f-p broj *)
fp[w_,k_]:=Module[{sw}, sw=MantissaExponent[w]; 
		                         sw[[1]]=Round[10^k sw[[1]] ]/10^k;
		                         sw[[1]] 10^sw[[2]]//N
		                                         ];

Primjer 3.1

a={{2,400},{1,1}}; b={200,1};
Print["a = ",MatrixForm[a],"   b = ",MatrixForm[b]]

(* Mathematica  - modul *)
Print["Mathematica - modul"]
x=LinearSolve[a,b]//N           (* Gauss *)

(* Gaussova metoda eliminacije *)
Print["Gaussova metoda eliminacije"]
Print["[A;b] = ",MatrixForm[A] ]
A={{2,400,200},{1,1,1}};
m=A[[2,1]]/A[[1,1]];
Do[
	A[[2,i]]=A[[2,i]]-m  A[[1,i]],
	{i,3}];
Print["A2 = ",MatrixForm[A]]
x[[2]]=A[[2,3]]/A[[2,2]];
x[[1]]=A[[1,3]]/A[[1,1]]-(A[[1,2]]/A[[1,1]]) x[[2]];
Print["x = ",N[x]]

(* Gaussova metoda eliminacije u k-znamekastoj floating point aritmetici *)
Print["Gaussova metoda eliminacije u f-p aritmetici"]
k=2;
A=N[{{2,400,200},{1,1,1}}];
Print["[A;b] = ",MatrixForm[A] ]
m=fp[A[[2,1]]/A[[1,1]],k]
Do[
	A[[2,i]]=fp[A[[2,i]]-fp[m  A[[1,i]],k],k],
	{i,3}];
Print["A2 = ",MatrixForm[A]]
x[[2]]=fp[A[[2,3]]/A[[2,2]],k];
x[[1]]=fp[A[[1,3]]/A[[1,1]],k]-fp[(fp[A[[1,2]]/A[[1,1]],k]) x[[2]],k];
Print["x = ",N[x]]

Gaussova metoda eliminacije

n=4;
a={{1,1,0,3},{2,1,-1,1},{3,-1,-1,2},{-1,2,3,-1}};b={4,1,-3,4};
Print[MatrixForm[a],MatrixForm[b]]
Do[
Do[
m=a[[i,l]]/a[[l,l]];
Do[
a[[i,k]]=a[[i,k]]-m  a[[l,k]],{k,l,n}];
   b[[i]]=b[[i]]-m  b[[l]],{i,l+1,n}];
Print[MatrixForm[a],MatrixForm[b] ],
{l,n-1}]
BS[n,a,b]

Primjer 3.2

n=2;a={{2,400},{1.,1}}; b={200,1}; db={1,0}; 
Print["A = ",MatrixForm[a],"    b = ",MatrixForm[b],"   db = ",MatrixForm[db]  ]
x=LinearSolve[a,b]//N; xtilda=LinearSolve[a,b+db];
Print["x = ",MatrixForm[x],"   xtilda = ",MatrixForm[xtilda]];
Delx=Max[Abs[xtilda-x]]; Delb=Max[Abs[db]];
Print["Delb = ",Delb,"    Delx = ",Delx,"    Delx/Delb = ", Delx/Delb];
Print["delb = ",Delb/Max[Abs[b]],"    delx = ",Delx/Max[Abs[x]],"    delx/delb = ", Max[Abs[b]]Delx/(Delb Max[Abs[x]])];
na=norma[a,n];na1=norma[Inverse[a],n];
Print["||A||=",na,"   ||A^(-1)||=",na1,"    cond(A)=",na na1] ;

Podaci (primjeri matrica)

(* Pozitivno definitna matrica kod LPNK *)
n=4;t={1,2,3,4,5};
a=Table[Apply[Plus,t^(i-1) t^(j-1)],{i,n},{j,n}];
b=Table[i,{i,n}]; db=Table[0,{i,n}]; db[[1]]=-.05;
Print["A = ",MatrixForm[a],"    b = ",MatrixForm[b],"   db = ",MatrixForm[db]];
na=norma[a,n]; na1=norma[Inverse[a],n];
Print["||A||=",na,"   ||A^(-1)||=",na1,"    cond(A)=",na na1] ;

(* Hilbertova matrica reda n *)
n=5; a=Table[1/(i+j-1),{i,n},{j,n}];
b=Table[i/100,{i,n}]; db=Table[.0001 i,{i,n}]; db[[1]]=-.005;
Print["A = ",MatrixForm[a],"    b = ",MatrixForm[b],"   db = ",MatrixForm[db]];
na=norma[a,n]; na1=norma[Inverse[a],n];
Print["||A||=",na,"   ||A^(-1)||=",na1,"    cond(A)=",na na1] ;

Cholesky dekompozicija pozitivcno definitne matrice

Cholesky[a,n]
FS[n,L,b]; BS[n,Transpose[L],x];x

LU dekompozicija

LU[a,n]
FS[n,L,b];BS[n,U,x];x//N

Invertiranje matrice primjenom LU dekompozicije

n=3;
a={{1,4,7},{2,5,8},{3,6,11}};
Print["A=",MatrixForm[a]]
LU[a,n];
inv=Table[0,{i,n},{j,n}];
Do[FS[n,L,IdentityMatrix[n][[i]]];
	     BS[n,U,x]; inv[[i]]=x,{i,n}]
Print["A^{-1}=",MatrixForm[Transpose[inv]]]

Iteratvne metode

(* Jacobijeva metoda *)
Jacobi[a_,b_,it_]:=Module[{n=Length[b],i,k,xs,xn}, xs=Table[0,{i,n}];xn=xs;
Do[
        Do[
	             xn[[i]]=(b[[i]]-Apply[Plus,ReplacePart[a[[i]],0,{i}] xs])/a[[i,i]],
				      {i,n}];
       Print["It_",k," = ",N[xn]]; xs=xn,
	                  {k,it}]
                                                                                ]

(* Gauss-Seidelova metoda *)
GS[a_,b_,eps_]:=Module[{k=0,xs,xn,n=Length[b],it=25},xs=xn=Table[1,{i,n}];
While[
Do[
xn[[i]]=(b[[i]]-Sum[a[[i,j]]xn[[j]],{j,i-1}]-Sum[a[[i,j]]xs[[j]],{j,i+1,n}])/a[[i,i]],{i,n}];
Max[Abs[xn-xs]]>eps &&k <it, xs=xn;k=k+1];
Print["it_",k," = ",N[xn]]]

Primjer

n=3;a={{1,.25,.5},{.25,1,.25},{.2,.25,1}};
             b={-1,0,1}; xx=LinearSolve[a,b];
Print["A = ",MatrixForm[a],"  b = ",MatrixForm[b], ", x = ",MatrixForm[xx]]
Print["Jacobijeva metoda"];Jacobi[a,b,8];Print[""]
Print["Gauss-Seidelova metoda"];GS[a,b,.005]