WICHTIG: Der Betrieb von goMatlab.de wird privat finanziert fortgesetzt. - Mehr Infos...

Mein MATLAB Forum - goMatlab.de

Mein MATLAB Forum

 
Gast > Registrieren       Autologin?   

Partner:




Forum
      Option
[Erweitert]
  • Diese Seite per Mail weiterempfehlen
     


Gehe zu:  
Neues Thema eröffnen Neue Antwort erstellen

Brauche Hilfe zu GMRES

 

wsm
Forum-Newbie

Forum-Newbie


Beiträge: 2
Anmeldedatum: 09.05.09
Wohnort: ---
Version: ---
     Beitrag Verfasst am: 28.12.2009, 16:53     Titel: Brauche Hilfe zu GMRES
  Antworten mit Zitat      
Hallo,

ich habe gerade GMRES nach dem Buch von Hanke Seite 137 ff. programmiert.

Die Funktionen für Arnoldi_Hessenberg und Givens_Rotation funktionieren und sind an bekannten Beispielen getestet.

Ich habe wohl noch ein Verständnisproblem bei dem gesamten Algorithmus.

Deswegen brauche ich ein wenig Hilfe.

Ich würde mcih freuen, wenn jemand der das Verfahren Kennt mal drüber schaut.

Vielen Dank
Guten Rutsch

W.

Code:

clc
clear

b=[6; 15; 9];
d=[0;0;0];
x0=[0;0;0];

A=[1 2 3; 4 5 6; 1 5 3];
[m,n]=size(A);
A

for k=1: 5      %Welches k ist sinnvoll?
    k
    [H,beta,V_M]=arnoldi_hessenberg(b,x0,A);    %beta = norm(rk)
   
    d(1)=beta;
    for i=2:n+1
        d(i)=0;
    end
    d;
   [R,Q,d]=Givens_Rotation(H,d);
    R,Q;
   
    c1=b(1:n); %Vektoranpaasung von Vektor b
    z=inv(R)*c1; %Lösung von R*z=c1
    x0=x0+V_M*z %neues x bestimmen
end

function [R,Q,b] = Givens_Rotation(H,b)

'Beginn Givens'
[m,n]=size(H); %m=max Zeilenzahl, n=max Spaltenzahl
H0=H;
for j=1:n
    for i=j+1:m
        G=eye(m);
        H
        % c und s bestimmen        
        hjj=H(j,j);
        hij=H(i,j);
        if(abs(hjj)>=abs(hij))
            t=hij/abs(hjj);
            n0=sqrt(1+abs(t)^2);
            c=sign(hjj)/n0;
            s=t/n0;
        else
            t=hjj/abs(hij);
            n0=sqrt(1+abs(t)^2);
            c=t/n0;
            s=sign(hij)/n0;
        end
        G(i,i)=c; %Besetzen der Givens-Rotationsmatrix
        G(i,j)=-s;
        G(j,i)=s;
        G(j,j)=c;
        G;
        b=G*b;
        H=G*H;
    end
   
    for i=1:n %Matrix R anpassen zu n x n Matrix aus H
        for j=1:n
        R(i,j)=H(i,j);    
        end
    end
    R;
    H;
    Qq=R/H0;   %Entspricht dem Produkt der Gn*Gn-1*...*G1 = Q'
     for i=1:n %Matrix Q anpassen zu nxn Matrix
        for j=1:m
        Q(i,j)=Qq(i,j);    
        end
     end
'Ende Givens_Rotation'
end

%Arnoldi-Verfahren
function [H,beta,V_M]=arnoldi_hessenberg(b, x0, A)
Info='Beginn Arnoldi'
[m,n]=size(A);
if m ~= n %prüft, ob Matrix quadratisch ist
    error ('Die Matrix ist nicht quadratisch!')
   
elseif length(b) ~= n %prüft, ob Lösungsvektor b die richtige Größe hat
    error ('b besitzt nicht die richtige Größe!')
   
elseif length(x0) ~= n %prüft, ob Startvektor x0 die richtige Größe hat
    error ('x0 besitzt nicht die richtige Größe!')
end

  r0=b-A*x0;              %Berechnet Residuum
    beta=norm(r0);        %Berechnet Norm des Residuums
    v=r0/beta;            %Berechnet ersten Vektor v1=r0/norm(r0)
     
    for j=1:n             %j ist Spaltenindex, n ist maximale Spaltenzahl
        for i=1:m
            V_M(i,j)=v(i); %Vektor v wird in Matrix V_M (Spalte j) übertragen
        end
        P=V_M*V_M';       %Berechnung des Orthogonalprojektors P
        q=(eye(n)-P)*A*v; %Berchnung von q zur Berechnung des nächsten v (16.3)
        H(j+1,j)=norm(q); %(**) Besetzung von H für i=j+1 für spätere Hessenbergmatrix
        v=q/norm(q);      %neues v wird berechnet        
    end    
     
 
    for j=1:n
        for i=1:m
            if i<=j
                for i2=1:n
                    u1(i2)=V_M(i2,i);
                end
                for i2=1:n
                    u2(i2)=V_M(i2,j);
                end

      %Berechnung des oberen Dreiecks der Hessenbergmatrix
            H(i,j)=u1*A*u2';
            end
            %Untere Diagonale (i=j+1) wurde bereits direkt nach der
            %Berechnung von q besetzt. (**)
            if i>j+1
                H(i,j)=0;
            end
        end
    end
      V_M
      H
      beta
Info='Ende Arnoldi'
end

 
Private Nachricht senden Benutzer-Profile anzeigen


Gast



Beiträge: ---
Anmeldedatum: ---
Wohnort: ---
Version: ---
     Beitrag Verfasst am: 11.03.2014, 17:20     Titel: Kann jemand helfen?
  Antworten mit Zitat      
Wir haben gerade genau das selbe Problem. Kann jemand weiterhelfen?
 
Jan S
Moderator

Moderator


Beiträge: 11.057
Anmeldedatum: 08.07.10
Wohnort: Heidelberg
Version: 2009a, 2016b
     Beitrag Verfasst am: 17.03.2014, 21:31     Titel: Re: Kann jemand helfen?
  Antworten mit Zitat      
Hallo,

Was bedeutet "genau das selbe Problem"? Die Orginal-Mitteilung enthielt ja keine genau Frage:
Zitat:

Ich habe wohl noch ein Verständnisproblem bei dem gesamten Algorithmus.
Deswegen brauche ich ein wenig Hilfe.
Ich würde mcih freuen, wenn jemand der das Verfahren Kennt mal drüber schaut.

Eine Antwort ist deshalb noch nicht möglich.

Gruß, Jan
Private Nachricht senden Benutzer-Profile anzeigen
 
Neues Thema eröffnen Neue Antwort erstellen



Einstellungen und Berechtigungen
Beiträge der letzten Zeit anzeigen:

Du kannst Beiträge in dieses Forum schreiben.
Du kannst auf Beiträge in diesem Forum antworten.
Du kannst deine Beiträge in diesem Forum nicht bearbeiten.
Du kannst deine Beiträge in diesem Forum nicht löschen.
Du kannst an Umfragen in diesem Forum nicht mitmachen.
Du kannst Dateien in diesem Forum posten
Du kannst Dateien in diesem Forum herunterladen
.





 Impressum  | Nutzungsbedingungen  | Datenschutz | FAQ | goMatlab RSS Button RSS

Hosted by:


Copyright © 2007 - 2024 goMatlab.de | Dies ist keine offizielle Website der Firma The Mathworks

MATLAB, Simulink, Stateflow, Handle Graphics, Real-Time Workshop, SimBiology, SimHydraulics, SimEvents, and xPC TargetBox are registered trademarks and The MathWorks, the L-shaped membrane logo, and Embedded MATLAB are trademarks of The MathWorks, Inc.