Matlab code pcg

Moderators: jkien, Xilvo

Reageer
Berichten: 4.246

Matlab code pcg

[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?
Quitters never win and winners never quit.

Berichten: 7.068

Re: Matlab code pcg

Ik maak hier ergens een fout:
Hoe kom je tot deze uitspraak? Krijg je een foutmelding? Doet de code niet wat je verwacht?

Berichten: 4.246

Re: Matlab code pcg

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.
Quitters never win and winners never quit.

Berichten: 7.068

Re: Matlab code pcg

Wat is M? Wat is b?

Berichten: 4.246

Re: Matlab code pcg

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))).
Quitters never win and winners never quit.

Berichten: 7.068

Re: Matlab code pcg

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 :D ).

Gebruikersavatar
Berichten: 7.224

Re: Matlab code pcg

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.
If I have seen further it is by standing on the shoulders of giants.-- Isaac Newton

Berichten: 4.246

Re: Matlab code pcg

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 :D

@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
Quitters never win and winners never quit.

Reageer