Pagina 1 van 1
Matlab code pcg
Geplaatst: za 27 dec 2008, 18:01
door dirkwb
[attachment=2957:1.PNG]
Ik maak hier ergens een fout:
Code: Selecteer alles
u = [ 0;zeros(n-2,1); 0]; % Let op: bij (preconditioned)CG moet deze nul zijn!
k = 0;
r = b;
while norm(r)>1e-12
z=inv(M)*r;
k = k+1
if k==1
p = z;
r_1 = b;
z_1 = z;
else
beta = (r_1'*z_1)/(r_2'*z_2);
p = z_1+ beta*p;
end
alpha = (r_1'*z_1)/(p'*A*p);
u = u + alpha*p;
r= r_1-alpha*A*p;
r_2 = r_1;
r_1 = r;
z_2 = z_1;
z_1 = z;
end
maar ik zie niet waar, kan iemand me helpen?
Re: Matlab code pcg
Geplaatst: ma 29 dec 2008, 00:42
door EvilBro
Ik maak hier ergens een fout:
Hoe kom je tot deze uitspraak? Krijg je een foutmelding? Doet de code niet wat je verwacht?
Re: Matlab code pcg
Geplaatst: zo 04 jan 2009, 22:36
door dirkwb
Krijg je een foutmelding? Doet de code niet wat je verwacht?
Nee, ja.
De code is een algoritme dat moet convergeren maar dat doet hij niet.
Re: Matlab code pcg
Geplaatst: wo 07 jan 2009, 22:37
door EvilBro
Wat is M? Wat is b?
Re: Matlab code pcg
Geplaatst: wo 07 jan 2009, 22:44
door dirkwb
Het algoritme lost op: Au=b maar om dit proces te versnellen bij hele grote matrices wordt er gebruikgemaakt van een 'preconditioner'. In dit geval M=PPT met P =sqrt(diag(A))).
Re: Matlab code pcg
Geplaatst: do 08 jan 2009, 08:21
door EvilBro
Mijn vraag was niet duidelijk. Ik wilde graag weten welke waarden jij gebruikt, zodat ik je programma kan runnen om te zien wat het doet (dan heb ik A dus ook nodig
).
Re: Matlab code pcg
Geplaatst: do 08 jan 2009, 19:59
door Bart
kun je overigens schrijven als
Je lost immers een set van lineaire vergelijkingen op en dat doe je te nimmer met een matrix inverse. Door het op deze manier te schrijven pakt Matlab de beste oplosmethode en deze is aanzienlijk sneller dan een matrix inverse.
Re: Matlab code pcg
Geplaatst: vr 09 jan 2009, 11:52
door dirkwb
Je lost immers een set van lineaire vergelijkingen op en dat doe je te nimmer met een matrix inverse. Door het op deze manier te schrijven pakt Matlab de beste oplosmethode en deze is aanzienlijk sneller dan een matrix inverse.
Ja maar dat is hier niet erg aangezien het om een diagonaalmatrix gaat
@Evilbro:
Code: Selecteer alles
%initialisering
A = zeros(n);
for j=1:n
A(j,j) = 2;
end
for j=1:n-1
A(j,j+1) = -1;
A(j+1,j) = -1;
end
h=1/(n+1);
A=A*1/h^2;% vergeet niet door h^2 te delen.
b=[ 1e3*ones(1,n/2) ones(1,n/2)]';%let op 1e3 en niet 1e-3
