Megjelenés: 2011
Oldalszám: 421 oldal
Formátum: B/5
ISBN: 978-963-2794-40-2
Témakör: Informatika, Matematika felsőfokon
Sorozat: Alkalmazott matematika
Eredeti ár: 4200 Ft
Webshop ár: 3150 Ft
KOSÁRBA
Programok
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% inkoga %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{verbatim}
function y=inkogaus(a,b,idd1,idd2,n,maxit,epsi,S)
% Kiszamitja az inkomplett Gauss-eliminaciot, majd iteracio-
% val oldja meg az Sy=b rendszert. Itt S az a matrix, amely
% a spdiags(a,idd1,n,n) utasitassal keszithet"o. Ennek meg-
% felel"oen a idd1 megadja a foglalt atlok helyet.
% Hivasa:
% y=inkogaus(a,b,idd1,idd2,n,maxit,epsi,S)
% ahol:
% y a lin. rendszer k"ozelit"o megoldasa,
% a az atlonkent tarolt matrix,
% b a jobboldali (oszlop-)vektor,
% idd1 az atlok helyet megado (oszlop-)vektor,
% idd2 az iteracio jobb konvergencia celjabol el"oirhato,
% kiegeszit"o atlok helye, mint oszlopvektor
% n az a matrix es a b vektor merete
% maxit a maximalis iteracioszam
% epsi a kilepesi feltetel abszolut pontossagi korlatja
% S kihagyhato: S=spdiags(a,idd1,n,n)
% Ha ez a matrix meg nincs, akkor inkoga kesziti.
[dda dda1]=size(idd1);
[idd,indd]=sort([idd1' idd2']); % a teljes indexvektor
%idd,indd,dda % "osszeallitasa
[dda1 dd]=size(idd);idi=0;
for i=1:dd
if idd(i)==0 % a f"oatlo helyenek
idi=i; % meghatarozasa
break
end
end
%idi
if idi==0
inputhiba=idi
else
idip=idi+1;idim=idi-1;
ddd=dd-idi;
for j=1:idim
k=j+1;llij=idd(j);
for i=1:ddd
iae=-1;iai=idd(i+idi)+llij;
for l=k:dd % az itt kesz"ul"o ia
if idd(l)==iai % indext"omb megadja azon
k=l+1; % atlok indexet, amelyen
iae=1; % valtozhatnak a matrix
ia(j,i)=l; % elemei (a f"oatlo alatti
break % atlokrol nezve)
elseif idd(l)>iai
break
end
end
if iae==-1
ia(j,i)=0;
end
end
end
% ia % elkesz"ult az indext"omb
id1=1; % bef"uzve a kiegeszit"o
for id=1:idi % nulla oszlopokat, olyan
if indd(id)>dda % alakra hozzuk a matrixot,
aa(:,id)=zeros(n,1); % hogy a felbontasnal keves
else % i>n-fele teszt kelljen
aa(:,id)=[zeros(-idd(id),1);a(1:n+idd(id),id1);];
id1=id1+1;
end
end
for id=idip:dd
if indd(id)>dda
aa(:,id)=zeros(n,1);
else
aa(:,id)=[a(1+idd(id):n,id1);zeros(idd(id),1)];
id1=id1+1;
end
end
clear indd
%aa % ezen matrix helyen
% t"ortenik a felbontas
zerus=0;
for i=1:n
if aa(i,idi)==0
zerus=i
break
else
d(i)=1/aa(i,idi);
for j=1:idim
ij=i-idd(j);
if ij<=n
llij=aa(ij,j)*d(i);
ll(ij,j)=llij;
for k=1:ddd % a felbontas legbels"obb
l=ia(j,k); % ciklusa
if l>0
aa(ij,l)=aa(ij,l)-llij*aa(i,idi+k);
end
end
end
end % j-ciklus vege
for id=1:ddd
uu(i,id)=aa(i,id+idi);
end
end
end
%aa,d,uu,ll % a felbontas matrixai
clear aa ia % elkesz"ultek
if zerus==0
y=b;
for i=1:n-1 % el"ore a nulladik
for j=1:idim % approximaciohoz
ij=i-idd(j);
if ij<=n
y(ij)=y(ij)-ll(ij,j)*y(i);
end
end
end
for i=n:-1:1 % vissza
w=y(i);
for id=idip:dd
is=i+idd(id);
if is<=n
w=w-uu(i,id-idi)*y(is);
end
end
y(i)=w*d(i);
end
if nargin<8
S=spdiags(a,idd1,n,n);
end
for it=1:maxit % az inkomplett felbontas
% utani iteracio kezdete, a
f=b-S*y; % maradekvektor kiszamitasa
z=norm(f,inf);
if z<epsi
break
end
for i=1:n-1 % el"ore
for j=1:idim
ij=i-idd(j);
if ij<=n
f(ij)=f(ij)-ll(ij,j)*f(i);
end
end
end
for i=n:-1:1 % vissza
w=f(i);
for id=idip:dd
is=i+idd(id);
if is<=n
w=w-uu(i,id-idi)*f(is);
end
end
f(i)=w*d(i);
end
y=y+f;
end % iteracio vege
norm_y=z,it
end % nullaoszto if vege
end % inputhiba if vege
\end{verbatim}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% qfigset %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{verbatim}
function qfigset(width,height,xlabels,subplab,linew)
%
%ke'szitette: Kolla'r Istva'n (BME) e's Gergo' Lajos (ELTE)
%
%
%a qfigset a rajzot (width x height) mereture allitja es
% beallitja a vonalvastagsagot linew ertekre
% width = az abra szelessege, alapertelmezes: az
% aranymetszes szerint a magassag ertekehez
% illesztve
% height = az abra magassaga, az aranymetszes szabalya
% szerint a szelesseg alapjan szamitva vagy
% 3.5 inch, ha a szelesseg nincs megadva
% xlabels = a keret es az x-tengely kozotti teruletet
% szabalyozza
% 'xtick': megengedi, hogy az x-tengely
% beosztasait megcimkezzuk
% 'nolabels': nem enged cimkezest
% 'plot' (alapertelmezes): megengedi az x-tengely
% es a beosztas megcimkezeset is
% 'subplot': megengedi az x-tengely es a beosztas
% megcimkezeset
% subplab = 'felirat' ha elozoleg a 'subplot' parameter
% volt megadva, akkor az x-tengely
% kozepere helyezi a feliratot.
% A keret meretenek az alapertelmezese 3.5 x 5.66 inch.
%
if nargin<5, linew=1.7; end
if nargin<3, xlabels='plot'; end
if strcmp(xlabels,'nolabels'), fcdist=0.1;
elseif strcmp(xlabels,'xtick'), fcdist=0.25;
elseif strcmp(xlabels,'plot'), fcdist=0.48;
elseif strcmp(xlabels,'subplot'), fcdist=0.70;
else error(['xlabels = ''',xlabels,''' is not valid'])
end
if nargin<4, subplab=''; end
if nargin<1, width = NaN; end
if nargin<2, height = NaN; end
gold=(sqrt(5)-1)/2;
if isempty(height), height=NaN; end
if isempty(width), width=NaN; end
%
%
% ha nem adtak meg szelesseget vagy magassagot
% (vagy egyiket sem),akkor beallitjuk az
% aranymetszes szerint (3.5 inch magassaggal)
%
%
if isnan(height) & isnan(width), height=3.5; end
if isnan(height), height=gold*width; end
if isnan(width), width=height/gold; end
%
Hp = 11; Wp = 8.4;
eval('tp=get(0,''TerminalProtocol'');','tp=''PCWIN'';')
if strcmp(tp,'PCWIN')
zunits=get(0,'Units');
set(0,'Units', 'inches')
scrdim=get(0,'ScreenSize');
if ~all(finite(scrdim))
scrdim
figure(1), scrdim=get(0,'ScreenSize')
end
set(0,'Units',zunits)
elseif ~strcmp(tp,'none') %ha csak rajz keszulhet
zunits=get(0,'Units');
set(0,'Units', 'inches')
scrdim=get(0,'ScreenSize');
set(0,'Units',zunits)
else
scrdim=get(0,'ScreenSize');
end
%
%
% itt szamitjuk ki a kepernyo meretet es a papir
% meretet a papir legfeljebb 8.4x11 inch meretu
% lehet (kb A4-es meret)
%
Hm=scrdim(4); Wm=scrdim(3);
Hp=min(Hp,Hm); Wp=min(Wm,Wp);
TopMarg = 0.1; BotMarg = fcdist;
LefMarg = 0.25; RiMarg = 0.25;
%
% ha vannak mar grafikus ablakok a kepernyon, es
% az aktualis ablak 'NextPlot' parametere 'new',
% akkor uj ablakot nyitunk
%
hc=get(0,'Children');
if ~isempty(hc)
if strcmp(get(gcf,'NextPlot'),'new'),
figure(max(hc)+1);
end
end
funits=get(gcf,'Units');
set( gcf, 'Units', 'inches' )
set( gcf, 'PaperUnits', 'inches' )
pwidth = width+RiMarg+LefMarg+0.5;
pheight = height+TopMarg+BotMarg+0.35;
% a grafikus ablak meretenek a beallitasa a
% kepernyon, az ablakot igyekszunk a jobb
% felso sarokba tenni
%
fleft = Wm - pwidth-0.05;
fbottom = Hm - pheight-0.7;
if fleft<0, pwidth=pwidth-fleft; fleft=0; end
if pwidth>Wp, pwidth=Wp; end
if pheight>Hp, pheight=Hp; end
if fbottom<0,pheight=pheight-fbottom;fbottom=0;end
set( gcf, 'Position',...
[ fleft fbottom pwidth pheight ])
% az abra meretenek a beallitasa a papiron
left = (Wp-pwidth);
bottom = Hp-pheight;
set( gcf, 'PaperPosition',...
[left-1.0 bottom pwidth pheight])
set(gcf,'DefaultAxesLineWidth',0.7)
set(gca,'LineWidth',0.7)
set(gcf,'DefaultLineLineWidth',linew)
set(gca,'DefaultLineLineWidth',linew)
%
% a tengelyek es a keret mereteinek a beallitasa
%
aunits=get(gca,'Units');
set( gca, 'Units', 'inches')
set( gca, 'Position', [ 0.75 fcdist width height ])
set(gcf,'DefaultAxesUnits','inches')
set( gcf, 'Units',funits)
set(gcf,'DefaultAxesPosition',[ 0.75/pwidth ...
fcdist/pheight width/pwidth height/pheight ])
set( gca, 'Units',aunits)
figure(gcf)
%
%
% a szoveg elhelyezese az x-tengely ala
%
%
if strcmp(xlabels,'subplot')
if ~isempty(subplab)
ax=axis;
xc=mean(ax(1:2));
yc=ax(3)-fcdist/pheight;
delete(findobj(gca,'Type','text','string',subplab))
text(xc,yc,subplab,...
'Horizontalalignment',...
'center','Verticalalignment','top');
end
end
%vege a qfigset programnak
\end{verbatim}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% gaufile %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{verbatim}
function [x,res]=gaufile(AF,AM,b,me,ep,imax)
%Nagymeretu, n ismeretlenes Ax=b egyenletrendszer megoldasa
%LU-modszerrel, iterativ javitassal
% Parameterek:
% AF : az A egyutthatomatrix elemeit oszlopfolytonos
% sorrendben tartalmazo file-nev (string)
% AM : segedfile-nev, a permutalt A matrix LU-felbontasat
% fogja tartalmazni a vegen (string)
% b : az egyenletrendszer jobboldala
% me : hany egyutthato fer be egyszerre a memoriaba,
% me>=n (integer)
% ep : a relativ modositas tureshatara (double)
%imax : az iteraciok maximalis szama (integer)
% x : az egyenletrendszer megoldasa
% res : a rezidualis hiba (res=Ax-b)
n=length(b);
% A memoriaban egyszerre tarthato oszlopok szama: kc
kc=floor(me/n);
if kc>n, kc=n;
elseif kc<1,
disp('HIBA! Kotelezoen: me>=egyenletek szama'), return
end
% Az A egyutthatomatrix particioinak szama: mc
mc=floor(n/kc);
if mc*kc>n, mc=mc-1;
elseif kc*mc<n, mc=mc+1; end
x=zeros(n,1); y=zeros(n,1);
% Az egyenletek permutaciojat mutato pointer vektor: p
p=1:n;
am=fopen(AM,'wb'); fclose(am);
[af,hibajel]=fopen(AF,'rb');
if af<0 disp('adatfile nyitasi HIBA az LU-felbontas elott!')
uzenet=hibajel, return
end
l=0;
n1=kc;
jel=0;
jels=0;
%Az LU-felbontas mc reszletben:
for k1=1:mc
if k1==mc, n1=n-(mc-1)*kc; end
%Az egyutthatomatrix kovetkezo n1 db oszlopanak beolvasasa
%(n1=kc, ha nem az utolso particio, egyebkent a maradek):
jel=jel+8*n*n1;
A=fread(af,[n,n1],'double');
if jel~=ftell(af)
particio=k1, disp('adatolvasasi HIBA!')
uzenet=ferror(af), return
end
%%Elmaradt eliminaciok potlasa
%(szukseg van az elozo, mar kesz oszlopokra):
[am,hibajel]=fopen(AM,'rb');
if am<0
particio=k1, disp('segedfile nyitasi HIBA!')
uzenet=hibajel, return
end
sjel=0;
for k=1:l
sjel=sjel+8*n;
col=fread(am,n,'double');
if sjel~=ftell(am)
particio=k1, disp('segedolvasasi HIBA!')
uzenet=ferror(am), return
end
for i=k+1:n
r=col(p(i));
A(p(i),1:n1)=A(p(i),1:n1)-r*A(p(k),1:n1);
end;
end
fclose(am);
%A memoriaban levo oszlopok tovabbi eliminacioja:
for k=1:n1
[s,t]=max(abs(A(p(l+k:n),k)));
if abs(s)<1e-50
disp('Figyelem! Az egyutthatomatrix szingularis,...
vagy majdnem szingularis!')
t=l+t+k-1; if k<t, p([l+k,t])=p([t,l+k]); end;
for i=l+k+1:n
r=A(p(i),k)/A(p(l+k),k); A(p(i),k)=r;
if k<n1
A(p(i),k+1:n1)=A(p(i),k+1:n1)-r*A(p(l+k),k+1:n1);
end;
end;
end;
%Elkeszult oszlopok felvitele az AM file-ba:
[am,hibajel]=fopen(AM,'ab');
if am<0
particio=k1, disp('segedfile nyitasi HIBA!')
uzenet=hibajel, return
end
jels=jels+8*n*n1;
as=fwrite(am,A,'double');
if jels~=ftell(am)|as~=n*n1
particio=k1, disp('segedirasi hiba HIBA!')
uzenet=ferror(am), return
end
fclose(am);
l=l+kc;
end
fclose(af); clear A;
%Befejezodott a permutalt A matrix LU-felbontasa.
%Mind az L, mind az U az AM file-ban van.
%A p mutatja a sorok permutaciojat, b1 jelolje a permutalt b-t.
%Az LUx=b1 megoldasa az eredeti Ax=b megoldasaval egyezik.
%Oldjuk meg elobb az Ly=b1 egyenletet, majd az Ux=y egyenletet!
y=solvtril(AM,b,p,me);
x=solvtriu(AM,y,p,me);
%Az x javitasa iteracioval mindaddig amig a relativ modositas
%nagyobb epszilon-nal, de legfeljebb imax lepesben:
for k=1:imax
%Eloszor az res=Ax-b rezidualis hibat szamitjuk:
res=-b;
n1=kc;
l=0;
[af,hibajel]=fopen(AF,'rb');
if af<0 disp('adatfile nyitasi HIBA az iteracio elott!')
uzenet=hibajel, return
end
jel=0;
for k1=1:mc
if k1==mc, n1=n-(mc-1)*kc; end
jel=jel+8*n*n1;
A=fread(af,[n,n1],'double');
if jel~=ftell(af)
disp('adatolvasasi hiba HIBA az iteracioban!')
particio=k1, uzenet=ferror(af), return
end
res=res+A*x(l+1:l+n1);
l=l+kc;
end
fclose(af); clear A;
%Most az Az=res egyenletrendszert oldjuk meg,
%itt is felhasznalva a permutalt A matrix LU felbontasat:
y=solvtril(AM,res,p,me);
z=solvtriu(AM,y,p,me);
d=norm(z)/norm(x);
%A kozelito x modositasa:
x=x-z;
if d<=ep return, end
end
return
\end{verbatim}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% solvtril %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{verbatim}
function [x]=solvtril(AM,b,p,me)
%Lx=b1 megoldasa, ahol L az AM-ben oszlopfolytonosan tarolt
%matrix sorok szerinti permutaltjanak also haromszog resze,
%a foatloban 1-ek feltetelezve, tovabba
%b1 a b parameter megfelelo permutacioja.
%Parameterek:
% AM : file-nev (string)
% b : az Lx=b1 egyenletrendszer
% jobboldalanak inverz permutaltja
% p : a permutacio mutatoja (integer vektor)
% me : az U hany eleme lehet egyszerre a memoriaban
n=length(b);
x=zeros(n,1);
% A memoriaban egyszerre tarthato oszlopok szama: kc
kc=floor(me/n);
if kc>n, kc=n;
elseif kc<1,
disp('HIBA! Kotelezoen: me>=egyenletek szama'), return
end
% Az L egyutthatomatrix particioinak szama: mc
mc=floor(n/kc);
if mc*kc>n, mc=mc-1;
elseif kc*mc<n, mc=mc+1; end
l=0;
n1=kc;
[am,hibajel]=fopen(AM,'rb');
if am<0 disp('file nyitasi HIBA az also haromszogmatrixban!')
uzenet=hibajel, return
end
jel=0;
%Az egyenletrendszer megoldasa mc reszletben:
for k1=1:mc
if k1==mc
n1=n-(mc-1)*kc;
end
jel=jel+8*n*n1;
%Az L egyutthatomatrix kovetkezo n1 db oszlopanak beolvasasa
%(n1=kc, ha nem az utolso particio, egyebkent a maradek):
L=fread(am,[n,n1],'double'); %AM(1:n,l+1:l+n1);
if jel~=ftell(am)
particio=k1,
disp('adatolvasasi HIBA az also haromszogmatrixban!')
uzenet=ferror(am), return
end
for i=1:n1
x(l+i)=b(p(l+i));
b=b-x(l+i)*L(1:n,i);
end
clear L;
l=l+kc;
end
fclose(am);
return
\end{verbatim}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% solvtriu %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{verbatim}
function [x]=solvtriu(AM,b,p,me)
%Ux=b megoldasa, ahol U az AM-ben oszlopfolytonosan tarolt
%matrix sorok szerinti permutaltjanak felso haromszog resze
%Parameterek:
% AM : file-nev (string)
% b : az egyenletrendszer jobboldala
% p : a permutacio mutatoja (integer vektor)
% me : az U hany eleme lehet egyszerre a memoriaban (integer)
n=length(b);
x=zeros(n,1);
% A memoriaban egyszerre tarthato oszlopok szama: kc
kc=floor(me/n);
if kc>n, kc=n;
elseif kc<1,
disp('HIBA! Kotelezoen: me>=egyenletek szama'), return
end
% Az U egyutthatomatrix particioinak szama: mc
mc=floor(n/kc);
if mc*kc>n, mc=mc-1;
elseif kc*mc<n, mc=mc+1; end
l=n;
n1=kc;
[am,hibajel]=fopen(AM,'rb');
if am<0 disp('file nyitasi HIBA a felso haromszogmatrixban!')
uzenet=hibajel, return
end
jel=8*n*n;
%Az egyenletrendszer megoldasa mc reszletben:
for k1=1:mc
if k1==mc
n1=n-(mc-1)*kc;
end
fseek(am,8*n*(l-n1),-1);
%Az U egyutthatomatrix kovetkezo n1 db oszlopanak beolvasasa
%(n1=kc, ha nem az utolso particio, egyebkent a maradek):
U=fread(am,[n,n1],'double'); %AM(1:n,l-n1+1:l);
if jel~=ftell(am)
particio=k1,
disp('adatolvasasi HIBA a felso haromszog matrixban!')
uzenet=ferror(am), return
end
for i=1:n1
s=U(p(l-i+1),n1-i+1)
if abs(s)<1e-50
disp('Figyelem! Az egyutthatomatrix szingularis,...
vagy majdnem szingularis!')
x(l-i+1)=b(l-i+1)/s;
b=b-x(l-i+1)*U(p,n1-i+1);
end
clear U;
l=l-kc; jel=jel-8*n*n1;
end
fclose(am);
return
\end{verbatim}