Matlab code pcg

[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?

[EvilBro]

Ik maak hier ergens een fout:

Hoe kom je tot deze uitspraak? Krijg je een foutmelding? Doet de code niet wat je verwacht?

[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.

[EvilBro]

Wat is M? Wat is b?

[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=PP^T met P =sqrt(diag(A))).

[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 ).

[Bart]

Code: Selecteer alles

z=inv(M)*r;

kun je overigens schrijven als

Code: Selecteer alles

z=M \ r;

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.

[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