Springen naar inhoud

Matlab code pcg


  • Log in om te kunnen reageren

#1

dirkwb

    dirkwb


  • >1k berichten
  • 4172 berichten
  • Moderator

Geplaatst op 27 december 2008 - 18:01

[attachment=2957:1.PNG]


Ik maak hier ergens een fout:

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?

Veranderd door dirkwb, 27 december 2008 - 18:10

Quitters never win and winners never quit.

Dit forum kan gratis blijven vanwege banners als deze. Door te registeren zal de onderstaande banner overigens verdwijnen.

#2

EvilBro

    EvilBro


  • >5k berichten
  • 6703 berichten
  • VIP

Geplaatst op 29 december 2008 - 00:42

Ik maak hier ergens een fout:

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

#3

dirkwb

    dirkwb


  • >1k berichten
  • 4172 berichten
  • Moderator

Geplaatst op 04 januari 2009 - 22:36

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.

#4

EvilBro

    EvilBro


  • >5k berichten
  • 6703 berichten
  • VIP

Geplaatst op 07 januari 2009 - 22:37

Wat is M? Wat is b?

#5

dirkwb

    dirkwb


  • >1k berichten
  • 4172 berichten
  • Moderator

Geplaatst op 07 januari 2009 - 22:44

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.

#6

EvilBro

    EvilBro


  • >5k berichten
  • 6703 berichten
  • VIP

Geplaatst op 08 januari 2009 - 08:21

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

#7

Bart

    Bart


  • >5k berichten
  • 7224 berichten
  • VIP

Geplaatst op 08 januari 2009 - 19:59

z=inv(M)*r;

kun je overigens schrijven als

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

#8

dirkwb

    dirkwb


  • >1k berichten
  • 4172 berichten
  • Moderator

Geplaatst op 09 januari 2009 - 11:52

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

Veranderd door dirkwb, 09 januari 2009 - 11:52

Quitters never win and winners never quit.





0 gebruiker(s) lezen dit onderwerp

0 leden, 0 bezoekers, 0 anonieme gebruikers

Ook adverteren op onze website? Lees hier meer!

Gesponsorde vacatures

Vacatures