Lineares Gleichungssystem

Neue Frage »

Manuel20 Auf diesen Beitrag antworten »
Lineares Gleichungssystem
Guten Abend allerseits!

Ich hätte eine kleine Frage:
Wie könnte ein (Matlab-)Code aussehen, der das QR-Verfahren anwendet, um ein lineares Gleichungssystem zu lösen.
Dabei denke ich, dass die Eingabe aus der regulären Matrix und dem rechte-Seite Vektor besteht. ..und die Ausgabe sollte dann natürlich der Lösungsvektor sein.

Ich habe gesehen, dass auf dem Board auch eine Sammlung von solchen "Problemen" vorhanden ist.
Ich denke hier jetzt speziell an QRmitH.zip - ein super Programm, das ich ein wenig "ummodeln" wollte (dass man zB den Vektor b eingeben kann), was aber ein wenig schief lief...

Ich danke Euch ganz herzlich und wünsche einen schönen Abend!
tigerbine Auf diesen Beitrag antworten »
RE: Lineares Gleichungssystem
Wie der Name des Programms sagt, QR mit Householder. Das bezieht sich natürlich nur auf die Zerlegung der Matrix.

Wo lief das mit b denn schief?

Ich nehme mir mal ein Testbeipiel aus dem Board. Du müsstest mir nun sagen, wo du b einbauen willst. Das Programm liefert die Rechnungen, die man in Klausuren von Hand zeigen soll, um die QR-Zerlegung einer Matrix zu bestimmen. Man kann ja nun das Lösen des LGS noch hintendranhängen.

Welche Matrizen brauchen wir denn dann um x auszurechnen? Was wissen wir über die? Q ist ja wunderschön. Augenzwinkern Und das letzte LGS dann, wegen der Gestalt von R durch Rückwärtssubstitution.

code:
1:
2:
3:
4:
5:
6:
7:
8:
9:
10:
11:
12:
13:
14:
15:
16:
17:
18:
19:
20:
21:
22:
23:
24:
25:
26:
27:
28:
29:
30:
31:
32:
33:
34:
35:
36:
37:
38:
39:
40:
41:
42:
43:
44:
45:
46:
47:
48:
49:
50:
51:
52:
53:
54:
55:
56:
57:
58:
59:
60:
61:
62:
63:
64:
65:
66:
67:
68:
69:
70:
71:
72:
73:
74:
75:
76:
77:
78:
79:
80:
81:
82:
83:
84:
85:
86:
87:
88:
89:
90:
91:
92:
93:
94:
95:
96:
97:
98:
99:
100:
101:
102:
103:
104:
105:
106:
107:
108:
109:
110:
111:
112:
113:
114:
115:
116:
117:
118:
119:
120:
121:
122:
123:
124:
125:
126:
127:
128:
129:
130:
131:
132:
133:
134:
135:
136:
137:
138:
139:
140:
141:
142:
143:
144:
145:
146:
147:
148:
149:
150:
151:
152:
>> QRmitH
Es wird eine QR-Zerlegung mit Householder berechnet.
====================================================
---------------------------
| H=I-ß*uu^T              |
| u=x+sign(x1)*||x||_2*e1 |
| ß=2/(u^Tu)              |
---------------------------
 
Matrix A eingeben: [-3,2,0;4,-1,5;0,1,-1]

A0 =

    -3     2     0
     4    -1     5
     0     1    -1

 
Durchgang 1 
================

x =

    -3
     4
     0


signx =

    -1


e1 =

     1
     0
     0


u =

    -8
     4
     0


beta =

   0.02500000000000

 

H =

  -0.60000000000000   0.80000000000000                  0
   0.80000000000000   0.60000000000000                  0
                  0                  0   1.00000000000000


Qk =

  -0.60000000000000   0.80000000000000                  0
   0.80000000000000   0.60000000000000                  0
                  0                  0   1.00000000000000


Ak =

   5.00000000000000  -2.00000000000000   4.00000000000000
  -0.00000000000000   1.00000000000000   3.00000000000000
                  0   1.00000000000000  -1.00000000000000

 
Durchgang 2 
================

x =

     1
     1


signx =

     1


e1 =

     1
     0


u =

   2.41421356237309
   1.00000000000000


beta =

   0.29289321881345

 

H =

  -0.70710678118655  -0.70710678118655
  -0.70710678118655   0.70710678118655


Qk =

   1.00000000000000                  0                  0
                  0  -0.70710678118655  -0.70710678118655
                  0  -0.70710678118655   0.70710678118655


Ak =

   5.00000000000000  -2.00000000000000   4.00000000000000
   0.00000000000000  -1.41421356237309  -1.41421356237309
   0.00000000000000                  0  -2.82842712474619

 
 
Die QR-Zerlegung mit Householder
--------------------------------

A =

    -3     2     0
     4    -1     5
     0     1    -1


Q =

  -0.60000000000000  -0.56568542494924  -0.56568542494924
   0.80000000000000  -0.42426406871193  -0.42426406871193
                  0  -0.70710678118655   0.70710678118655


R =

   5.00000000000000  -2.00000000000000   4.00000000000000
   0.00000000000000  -1.41421356237309  -1.41421356237309
   0.00000000000000                  0  -2.82842712474619
 
 
Manuel20 Auf diesen Beitrag antworten »
RE: Lineares Gleichungssystem
Hallo tigerbine!

Ja, dass die QR-Zerlegung wunderbar klappt, habe ich gesehen - super! smile
Was ich eben noch 'einbauen' wollte, ist, dass man am Anfang nicht nur die Matrix A eingibt, sondern auch gleich der rechte-Seite-Vektor b. (alles aus dem komplexen Zahlenbereich)
Und am Schluss wollte ich, dass das Programm den Loesungsvektor x ausgibt..aber soweit hat der Code bei mir nicht funktioniert, weil ich wohl bei der Eingabe von b etwas falsch implementiert habe :S

..ich versuche mal, den Fehler zu finden. Generell die Frage: die Eingabe von b wuerdest duch auch direkt am Anfang, und nicht nach der vollstaendigen QR-Berechnung einfuegen, oder?
tigerbine Auf diesen Beitrag antworten »
RE: Lineares Gleichungssystem
Wann ist egal. Solange es geschieht, bevor du drauf zugreifst. Big Laugh

code:
1:
2:
3:
b=input('Vektor b eingeben: ');


Dann musst du eben nur überlegen, wie du aus QRx=b x berechnen kannst.
Manuel20 Auf diesen Beitrag antworten »
RE: Lineares Gleichungssystem
Das habe ich mir schon ueberlegt:

Man kann ja das Gleichungssystem QRx = Ax = b in zwei Stufen loesen:
1.) Man berechnet y = Q^Hb = Pb
2.) Man loest danach Rx = A^{n=1} *x = y auf.

Mir ist gerade etwas aufgefallen:
Vergiss das obere:
Es wuerde doch ganz simpel gehen, indem man QRx durch QR teilt, um das x zu bekommen, bzw. Ax durch A oder b durch A oder b durch QR smile
..wuerde das klappen? (bzw wie (codemaessig))
tigerbine Auf diesen Beitrag antworten »
RE: Lineares Gleichungssystem
Durch eine Matrix teilen.... aua.... Teufel

Die erste Idee ist die richtige. Und ich sagte das ja auch schon. Augenzwinkern

Zitat:
Q ist ja wunderschön. Augenzwinkern Und das letzte LGS dann, wegen der Gestalt von R durch Rückwärtssubstitution.
Manuel20 Auf diesen Beitrag antworten »
RE: Lineares Gleichungssystem
Ouw..das mit der Division durch eine Matrix hat *niemand* gesehen..
*peinoo*

Nun aber zum Code.
Den habe ich soweit (minim) verändert:

clear all;
disp('Es wird eine QR-Zerlegung mit Householder berechnet.')
disp('====================================================')
disp('---------------------------')
disp('| H=I-ß*uu^T |')
disp('| u=x+sign(x1)*||x||_2*e1 |')
disp('| ß=2/(u^Tu) |')
disp('---------------------------')
disp(' ')
A0=input('Matrix A eingeben: ')
b0=input('Vektor b eingeben: ');
disp(' ')
n=length(A0);
m=length(A0');
c=length(b);
A=A0;
b=b0;
min=min(m-1,n);
Q=eye(n);

for k=1:min
fprintf('Durchgang %g \n', k)
disp('================')

clear x;
clear u;
clear beta;
clear e1;

e1=zeros(m+1-k,1);
for i=k:min+1
x(i+1-k,1)=A(i,k);
end
x
x2=norm(x,2);
signx=sign(x(1,1))
e1(1,1)=1;
e1

u=x+signx.*x2.*e1
beta=2/(u'*u)

id=eye(n+1-k);
disp(' ')
H=id-beta.*(u*u')

Qk=eye(n);
for i=k:n
for j=k:n
Qk(i,j)=H(i+1-k,j+1-k);
end
end
Qk
Q=Qk*Q;
Ak=Qk*A
A=Ak;
disp(' ')
pause
end
R=A;
A=A0;
Q=Q';
x=Q/R;
disp(' ')
disp('Die QR-Zerlegung mit Householder')
disp('--------------------------------')
A
Q
R
x

Das Problem, das momentan noch besteht ist, wie ich den Vektor b als solchen definiere..
tigerbine Auf diesen Beitrag antworten »
RE: Lineares Gleichungssystem
Verstehe das Problem nicht.

Was soll das sein

code:
1:
2:
x=Q/R;


Ich sehe nirgends das doppelte Lösen eines LGS. Dann hättest du ja auch x=A\b eingeben können.
Manuel20 Auf diesen Beitrag antworten »
RE: Lineares Gleichungssystem
Zitat:
Original von tigerbine

Ich sehe nirgends das doppelte Lösen eines LGS. Dann hättest du ja auch x=A\b eingeben können.


Wie könnte denn das aussehen - codemässig?
tigerbine Auf diesen Beitrag antworten »
RE: Lineares Gleichungssystem
Ach komm, so was kannst du doch selbst nachschlagen. Ich sagte schon, das Q was besonderes ist und man hier wohl das tun wird, was man sonst nur tun wird: Invertieren

Ferner ist Rückwärtssubstitution ein dankbarer Suchbegriff, auch für die Boardsuche. Augenzwinkern

Das mit A\b war als Witz gemeint. Dann braucht man auch nicht doe QR-Zerlegung zu berechnen, wenn man den Matlablöser "\" nimmt. Augenzwinkern
Manuel20 Auf diesen Beitrag antworten »
RE: Lineares Gleichungssystem
Okey, ich habe nun folgendes gemacht:

code:
1:
2:
3:
4:
5:
6:
7:
8:
9:
10:
11:
12:
13:
14:
15:
16:
17:
18:
19:
20:
21:
22:
23:
24:
25:
26:
27:
28:
29:
30:
31:
32:
33:
34:
35:
36:
37:
38:
39:
40:
41:
42:
43:
44:
45:
46:
47:
48:
49:
50:
51:
52:
53:
54:
55:
56:
57:
58:
59:
60:
61:
62:
63:
64:
65:
66:
67:
68:
69:
70:
71:
72:
73:
74:
75:
76:
77:
78:
79:
80:
81:
82:
83:
84:
85:
86:
87:
88:
89:
90:
91:
92:
93:
94:
95:
96:
97:
98:
99:
%QR-Zerlegung mit Householder

clear all;
disp('Es wird eine QR-Zerlegung mit Householder berechnet.')
disp('====================================================')
disp('---------------------------')
disp('| H=I-ß*uu^T |')
disp('| u=x+sign(x1)*||x||_2*e1 |')
disp('| ß=2/(u^Tu) |')
disp('---------------------------')
disp(' ')
A0=input('Matrix A eingeben: ')
b=input('Vektor b eingeben: ');
disp(' ')
disp(' ')
n=length(A0);
m=length(A0');
A=A0;
min=min(m-1,n);
Q=eye(n);
% Ist b ein Spaltenvektor ?
[q,p] = size(b);
if q~=n
error('b hat die falsche Dimension!');
end


for k=1:min
fprintf('Durchgang %g n', k)
disp('================')

clear x;
clear u;
clear beta;
clear e1;

e1=zeros(m+1-k,1);
for i=k:min+1
x(i+1-k,1)=A(i,k);
end
x
x2=norm(x,2);
signx=sign(x(1,1))
e1(1,1)=1;
e1

u=x+signx.*x2.*e1
beta=2/(u'*u)

id=eye(n+1-k);
disp(' ')
H=id-beta.*(u*u')

Qk=eye(n);
for i=k:n
for j=k:n
Qk(i,j)=H(i+1-k,j+1-k);
end
end
Qk
Q=Qk*Q;
Ak=Qk*A
A=Ak;
x=Ab;
disp(' ')
pause

% Loesen des Gleichungssystems

% Schritt 2: Vorwaertssubstitution Ay=b
y = zeros(n,1);
for j=1:n
summe = 0;
for k=1:j-1
summe = summe + A(j,k) * x(k);
end;
x(j) = (b(j) - summe) / A(j,j);
end


% Schritt 2: Rueckwaertssubstitution Rx=z mit R=A'
x = zeros(n,1);
for j=n:-1:1
summe = 0;
for k=j+1:n
summe = summe + R(j,k) * x(k);
end;
x(j) = (y(j)- summe) / R(j,j);
end



disp('Loesungsvektor x=');
disp(x);
end


Leider funktioniert's aber immer noch nicht ganz, wie gewünscht.Zwar kann ich nun den Vektor b eingeben, und ich denke, im Hintergrund wird auch das Richtige gemacht, aber die Ausgabe ist noch nicht so, wie gewünscht..Könntest du dir den Code evtl. einmal anschauen und mir sagen, was du wo anders machen würdest?Besten Dank und einen schönen Abend!
tigerbine Auf diesen Beitrag antworten »
RE: Lineares Gleichungssystem
Du weißt wie das mit Programmierübungen ist. Läuft oder läuft nicht. Da gibt es keine Korrektur. Big Laugh

Was heißt läuft nicht wie gewünscht? Berechne mit A\b die korrekte Lösung, wähle ganzzahliges Beispiel. Dann rechne nach QR von Hand, nimm die ";" raus und schaue, ab wo nicht korrekt gerechnet wird.
Manuel20 Auf diesen Beitrag antworten »
RE: Lineares Gleichungssystem
Hehe smile
Ich würde sagen, dass alles vorhanden wäre..
Aber ich seh den (die) Fehler nicht..:S
Jetzt erhalte ich sogar eine Warnung (bei nicht zu komplizierter Matrix A):

Warning: Matrix is close to singular or badly scaled.
Results may be inaccurate. RCOND = 1.967192e-017.
tigerbine Auf diesen Beitrag antworten »
RE: Lineares Gleichungssystem
Gib mal die Matrix und das b. Und du musst mir schon genauer sagen, was nicht passt.
Manuel20 Auf diesen Beitrag antworten »
RE: Lineares Gleichungssystem
Ich habe einfach ein Beispiel durgelassen:
A=[1,2,3;4,5,6;7,8,9]
und b = [7;11;12]
tigerbine Auf diesen Beitrag antworten »
RE: Lineares Gleichungssystem
Gib mal nur die Matrix ein und dann rank(A). Schlechtes Beispiel, was du genommen hast. Augenzwinkern
Manuel20 Auf diesen Beitrag antworten »
RE: Lineares Gleichungssystem
Tatsächlich: Rang = 2

Besseres Beispiel:
A=[32,13,34;23,92,19;23,40,13]
b=[3;6;2]

Hier erscheint kein Fehler, sondern folgendes:

code:
1:
Es wird eine QR-Zerlegung mit Householder berechnet.====================================================---------------------------| H=I-ß*uu^T || u=x+sign(x1)*||x||_2*e1 || ß=2/(u^Tu) |--------------------------- Matrix A eingeben: [32,13,34;23,92,19;23,40,13]A0 = 32 13 34 23 92 19 23 40 13Vektor b eingeben: [3;6;2] Durchgang 1 ================x = 32 23 23signx = 1e1 = 1 0 0u = 77.6289 23.0000 23.0000beta = 2.8232e-004 H = -0.7013 -0.5041 -0.5041 -0.5041 0.8507 -0.1493 -0.5041 -0.1493 0.8507Qk = -0.7013 -0.5041 -0.5041 -0.5041 0.8507 -0.1493 -0.5041 -0.1493 0.8507Ak = -45.6289 -75.6537 -39.9746 -0.0000 65.7336 -2.9173 -0.0000 13.7336 -8.9173

Der gesuchte Lösungsvektor x, der eigentlich auch noch ausgegeben werden sollte, erscheint leider aber nicht...
Manuel20 Auf diesen Beitrag antworten »
RE: Lineares Gleichungssystem
Hier nochmals die Ausgabe - ein wenig "leserlicher":


Es wird eine QR-Zerlegung mit Householder berechnet.
====================================================
---------------------------
| H=I-ß*uu^T |
| u=x+sign(x1)*||x||_2*e1 |
| ß=2/(u^Tu) |
---------------------------

Matrix A eingeben: [32,13,34;23,92,19;23,40,13]

A0 =

32 13 34
23 92 19
23 40 13

Vektor b eingeben: [3;6;2]


Durchgang 1
================

x =

32
23
23


signx =

1


e1 =

1
0
0


u =

77.6289
23.0000
23.0000


beta =

2.8232e-004



H =

-0.7013 -0.5041 -0.5041
-0.5041 0.8507 -0.1493
-0.5041 -0.1493 0.8507


Qk =

-0.7013 -0.5041 -0.5041
-0.5041 0.8507 -0.1493
-0.5041 -0.1493 0.8507


Ak =

-45.6289 -75.6537 -39.9746
-0.0000 65.7336 -2.9173
-0.0000 13.7336 -8.9173
tigerbine Auf diesen Beitrag antworten »
RE: Lineares Gleichungssystem
Was soll ich damit?

Die QR zerlegung läuft korrekt. Dann eben

* y=Q'b

* Rx=y mit Rückwärtssubstitution. Code steht doch im Workshop.

Läuft dann auch alles korrekt.
tigerbine Auf diesen Beitrag antworten »
RE: Lineares Gleichungssystem
code:
1:
2:
3:
4:
5:
6:
7:
8:
9:
10:
11:
12:
13:
14:
15:
16:
17:
18:
19:
20:
21:
22:
23:
24:
25:
26:
% Loesen des Gleichungssystems

% Schritt 2: Vorwaertssubstitution Ay=b
y = zeros(n,1);
for j=1:n
summe = 0;
for k=1:j-1
summe = summe + A(j,k) * x(k);
end;
x(j) = (b(j) - summe) / A(j,j);
end


% Schritt 2: Rueckwaertssubstitution Rx=z mit R=A'
x = zeros(n,1);
for j=n:-1:1
summe = 0;
for k=j+1:n
summe = summe + R(j,k) * x(k);
end;
x(j) = (y(j)- summe) / R(j,j);
end





1. Warum A? Du sollst doch Q nehmen. Ferner geht Vorwärstsub doch nur mit einer unteren Dreiecksmatrix.

2. Die Schleife ist zu lang.

3. Du hast bei deinem VS schon x als Variable. Was soll das da? und bei der RS greifst du auf y zu, was du als Nullvektor programmiert hast. Und das x hast du auch wieder auf Null gesetzt.
Manuel20 Auf diesen Beitrag antworten »
RE: Lineares Gleichungssystem
Guten Tag!

So, ich habe nun das, was falsch war, rausgelöscht, und den Code vom Workshop implementiert:

code:
1:
2:
3:
4:
5:
6:
7:
8:
9:
10:
11:
12:
13:
14:
15:
16:
17:
18:
19:
20:
21:
22:
23:
24:
25:
26:
27:
28:
29:
30:
31:
32:
33:
34:
35:
36:
37:
38:
39:
40:
41:
42:
43:
44:
45:
46:
47:
48:
49:
50:
51:
52:
53:
54:
55:
56:
57:
58:
59:
60:
61:
62:
63:
64:
65:
66:
67:
68:
69:
70:
71:
72:
73:
74:
75:
76:
77:
78:
79:
80:
81:
82:
83:
84:
85:
86:
87:
88:
89:
90:
91:
92:
93:
94:
95:
96:
97:
98:
99:
100:
101:
102:
103:
104:
105:
106:
107:
108:
109:
110:
111:
%QR-Zerlegung mit Householder

clear all;
disp('Es wird eine QR-Zerlegung mit Householder berechnet.')
disp('====================================================')
disp('---------------------------')
disp('| H=I-ß*uu^T |')
disp('| u=x+sign(x1)*||x||_2*e1 |')
disp('| ß=2/(u^Tu) |')
disp('---------------------------')
disp(' ')
A0=input('Matrix A eingeben: ')
b=input('Vektor b eingeben: ');
disp(' ')
disp(' ')
n=length(A0);
m=length(A0');
A=A0;
min=min(m-1,n);
Q=eye(n);
% Ist b ein Spaltenvektor ?
[q,p] = size(b);
if q~=n
error('b hat die falsche Dimension!');
end


for k=1:min
fprintf('Durchgang %g n', k)
disp('================')

clear x;
clear u;
clear beta;
clear e1;

e1=zeros(m+1-k,1);
for i=k:min+1
x(i+1-k,1)=A(i,k);
end
x
x2=norm(x,2);
signx=sign(x(1,1))
e1(1,1)=1;
e1

u=x+signx.*x2.*e1
beta=2/(u'*u)

id=eye(n+1-k);
disp(' ')
H=id-beta.*(u*u')

Qk=eye(n);
for i=k:n
for j=k:n
Qk(i,j)=H(i+1-k,j+1-k);
end
end
Qk
Q=Qk*Q;
Ak=Qk*A
A=Ak;
disp(' ')
pause
end
% Loesen des Gleichungssystems

% Schritt 2: Vorwaertssubstitution
function f = forward(L,b)
[n,n]=size(L);
L
b
%Lx=b Lösung eines GLS mit unterer Dreiecksmatrix
x(1,1)=b(1)./L(1,1);
for i=2:n
sum=0;
for j=1:(i-1)
sum=sum+L(i,j).*x(j,1);
end
x(i,1)=(b(i)-sum)./(L(i,i));
end
x

%Probe
Lb


% Schritt 2: Rueckwaertssubstitution
function b = backward(R,u)
[n,n]=size(R);

%Rx=u Lösung eines GLS mit oberer Dreiecksmatrix
x(n,1)=u(n,1)./R(n,n);
for i=(n-1):-1:1
sum=0;
for j=(i+1):n
sum=sum+R(i,j).*x(j,1);
end
x(i,1)=(u(i)-sum)./(R(i,i));
end
x

%Probe
%Ru


disp('Loesungsvektor x=');
disp(x);


Meine Fragen: Ich habe ja ein lineares Gleichungssystem: Ax = b
Für das Programm sollte man die Matrix A und den rechte-Seite-Vektor b eingeben sollen.
--> A kann man eingeben, genau so wie b.
Gut. Am Schluss sollte das Programm den Lösungsvektor x ausgeben, was leider nicht geschieht.
Hierzu habe ich Fragen:
1.)Bei der Vorwärts- und Rückwärtssubstitution: Sind dort die Variablen richtig gewählt, also beziehen sie sich auf das, was eigentlich "gemeint" ist?
2.) Zum Beispiel bin ich nicht sicher, aber x wird doch eigentlich anders verwendet. Was ändere ich dann (am besten) wo, und zu was?
3.) Die Ausgabe von x (also dem Lösungsvektor) würde man doch so machen:
code:
1:
2:
3:
disp('Loesungsvektor x=');
disp(x);

oder nicht?
tigerbine Auf diesen Beitrag antworten »
RE: Lineares Gleichungssystem
In Deutsch hieße es Themaverfehlung. Wir berechnen eine QR-Zerlegung und nicht eine LR-Zerlegung. Ich sagte doch gestern schon, dass eine Vorwärtssubstitution hier fehl am Platz ist.

Lies doch bitte, was man für dich schreibt.

Zitat:

1.)Bei der Vorwärts- und Rückwärtssubstitution: Sind dort die Variablen richtig gewählt, also beziehen sie sich auf das, was eigentlich "gemeint" ist?


Es ist schon falsch, forward zu machen.

Zitat:

2.) Zum Beispiel bin ich nicht sicher, aber x wird doch eigentlich anders verwendet. Was ändere ich dann (am besten) wo, und zu was?


Verstehe die Frage nicht. x ist der Vektor, der am Ende rauskommt.

Zitat:

3.) Die Ausgabe von x (also dem Lösungsvektor) würde man doch so machen:


Egal. Da du im Code schon "x" (ohne Augenzwinkern stehen hast, gibt man x eh schon aus. Du kannst das optisch so aufblähen, wie du magst. Augenzwinkern


Ich wiederhole mich ein letztes mal. Das zu programmieren dauert max 5 Minuten. Wir werden das doch heute noch gebacken bekommen....

Zitat:

* y=Q'b

* Rx=y mit Rückwärtssubstitution. Code steht doch im Workshop.


Auf einen erfolgreichen nächsten post. Prost
Manuel20 Auf diesen Beitrag antworten »
RE: Lineares Gleichungssystem
Oke, dass die Vorwärtssubstitution hier fehl am Platz ist, habe ich mittlerweile eingesehen (endlich :P)

Wie meinst du das aber:
* y=Q'b

* Rx=y mit Rückwärtssubstitution.

Also wo soll ich das einfügen, und was soll Q' bedeuten?

Das wären noch meine zwei Fragen, damit der nächste Post dann erfolgreich wird Augenzwinkern
tigerbine Auf diesen Beitrag antworten »
RE: Lineares Gleichungssystem
Wenigstens kommt die Einsicht so langsam. Ich hab schon Fusseln am Mund.

Qy=b <=> y=Q'b

Warum? Benutzte, dass Q orthogonal. Dann weißt du was Q' ist und auch schon, wie man es in matlab eingibt.

Einfügen? Na alles nach meinem Algorithmus.
Manuel20 Auf diesen Beitrag antworten »
RE: Lineares Gleichungssystem
Zitat:
Original von tigerbine

Dann weißt du was Q' ist und auch schon, wie man es in matlab eingibt.



..nicht wirklich.. :S
Also y ist eine neue, Variable, oder?
Theoretisch wäre das aber das x respektive die Variable, die ich am Schluss ausgegeben haben will...(oder?)
Mann..im Moment habe ich ein echtes Variabel-Durcheinander.. traurig
tigerbine Auf diesen Beitrag antworten »
RE: Lineares Gleichungssystem
Ich versteh einfach nicht, was du nicht verstehst. Tränen

Wir wollen lösen Ax=b. Wir zerlegen A=QR. Das macht dann QRx=b. Nun lösen wir 2 Teilprobleme. Q(Rx)=b, also Qy=b mit Rx=y. Und Q ist verdammt nochmal orthogonal. Warum ist das so eine tolle Information? Weil man dann die Inverse so einfach berechnen kann. Das ist nämlich die ..... von Q. Und die heißt in matlab Q'.
Manuel20 Auf diesen Beitrag antworten »
RE: Lineares Gleichungssystem
ahhhh..jetzt ist vieles klarer, zumindest was das Verständnis angeht smile
Weil Q orthogonal ist, lassen sich die Inverse einfach berechnen, was dann die (komplex) transponierte Matrix der Matrix Q ist smile

Okey, dass diese in matlab mit Q' bezeichnet wird, wusst' ich eben nicht. (Ich "arbeite" erst seit ca. 2, 3 Wochen mit dem Programm..)

Super, das Problem ist nun vollständig verstanden.

Nur happert's nach wie vor mit der Implementierung in matlab.
Die QR-Zerlegung funktioniert ja einwandfrei.
Den Vektor b kann man ebenfalls eingeben. Da die QR-Zerlegung funktioniert, verändere ich also den Teil ab "for k=1:min" bis zu dessen end nicht (bis auf "pause", das hab ich entfernt).
Nach diesem end füge ich folgendes an:
code:
1:
% Loesen des Gleichungssystems% Rueckwaertssubstitution function b = backward(R,u)[n,n]=size(R);%Rx=u Lösung eines GLS mit oberer Dreiecksmatrixx(n,1)=u(n,1)./R(n,n); for i=(n-1):-1:1 sum=0; for j=(i+1):n sum=sum+R(i,j).*x(j,1); end x(i,1)=(u(i)-sum)./(R(i,i));endx%Probe%R\uRx=y;Qy=b;Q';disp(' ')

Es folgt am Schluss, nach Eingabe des Vektors b und der Matrix A die Fehlermeldung:
??? function b = backward(R,u)
|
Error: Function definitions are not permitted at the prompt or in scripts.
Manuel20 Auf diesen Beitrag antworten »
RE: Lineares Gleichungssystem
Sorry:
code:
1:
2:
3:
4:
5:
6:
7:
8:
9:
10:
11:
12:
13:
14:
15:
16:
17:
18:
19:
20:
21:
22:
23:
24:
% Loesen des Gleichungssystems
% Rueckwaertssubstitution
function b = backward(R,u)
[n,n]=size(R);

%Rx=u Lösung eines GLS mit oberer Dreiecksmatrix
x(n,1)=u(n,1)./R(n,n);
for i=(n-1):-1:1
sum=0;
for j=(i+1):n
sum=sum+R(i,j).*x(j,1);
end
x(i,1)=(u(i)-sum)./(R(i,i));
end
x

%Probe
%Ru
Rx=y;
Qy=b;
Q';
disp(' ')
tigerbine Auf diesen Beitrag antworten »
RE: Lineares Gleichungssystem
Ergänze das Programm um

code:
1:
2:
3:
4:
5:
6:
7:
8:
9:
10:
11:
12:
13:
14:
15:
16:
17:
18:
19:
20:
21:
22:
23:
24:
25:
26:
%Berechnung einer LGS Lösung
 disp(' ');
 ans=input('Soll noch ein LGS gelöst werden? 0-nein, 1-ja: ');
 
 if ans==0
     break
 else
 b=input('Vektor b eingeben: ')
 
 %Qy=b lösen. Q ist orthogonal Q^(-1)=Q'
   y=Q'*b;
 
 %Rx=y durch Rückwärtssubstitution lösen.
 x(n,1)=y(n,1)/R(n,n);
 sum=0;
 for i=(n-1):-1:1
     for j=i+1:n         
     sum=sum+R(i,j)*x(j,1);
     end
     x(i,1)=(y(i,1)-sum)/R(i,i);
     sum=0;
 end
 x
 
 %Kontrolle
 xm=A\b


functions benutzte ich bei solchen Programmen nicht oft.
Manuel20 Auf diesen Beitrag antworten »
RE: Lineares Gleichungssystem
Wow, du bist ein echtes Genie! smile

Also schau, der vollständige Code sähe nun wie folgt aus:

code:
1:
2:
3:
4:
5:
6:
7:
8:
9:
10:
11:
12:
13:
14:
15:
16:
17:
18:
19:
20:
21:
22:
23:
24:
25:
26:
27:
28:
29:
30:
31:
32:
33:
34:
35:
36:
37:
38:
39:
40:
41:
42:
43:
44:
45:
46:
47:
48:
49:
50:
51:
52:
53:
54:
55:
56:
57:
58:
59:
60:
61:
62:
63:
64:
65:
66:
67:
68:
69:
70:
71:
72:
73:
74:
75:
76:
77:
78:
79:
80:
81:
82:
83:
84:
85:
86:
87:
88:
89:
90:
91:
92:
93:
94:
95:
96:
97:
98:
%QR-Zerlegung mit Householder

clear all;
disp('Es wird eine QR-Zerlegung mit Householder berechnet.')
disp('====================================================')
disp('---------------------------')
disp('| H=I-ß*uu^T |')
disp('| u=x+sign(x1)*||x||_2*e1 |')
disp('| ß=2/(u^Tu) |')
disp('---------------------------')
disp(' ')
A0=input('Matrix A eingeben: ')
b=input('Vektor b eingeben: ');
disp(' ')
disp(' ')
n=length(A0);
m=length(A0');
A=A0;
min=min(m-1,n);
Q=eye(n);

% Ist b ein Spaltenvektor ?
[q,p] = size(b);
if q~=n
error('b hat die falsche Dimension!');
end


for k=1:min
fprintf('Durchgang %g n', k)
disp('================')

clear x;
clear u;
clear beta;
clear e1;

e1=zeros(m+1-k,1);
for i=k:min+1
x(i+1-k,1)=A(i,k);
end
x
x2=norm(x,2);
signx=sign(x(1,1))
e1(1,1)=1;
e1

u=x+signx.*x2.*e1
beta=2/(u'*u)

id=eye(n+1-k);
disp(' ')
H=id-beta.*(u*u')

Qk=eye(n);
for i=k:n
for j=k:n
Qk(i,j)=H(i+1-k,j+1-k);
end
end
Qk
Q=Qk*Q;
Ak=Qk*A
A=Ak;
disp(' ')


end


%Berechnung einer LGS Lösung
disp(' ');
ans=input('Soll noch ein LGS gelöst werden? 0-nein, 1-ja: ');

if ans==0
break
else
b=input('Vektor b eingeben: ')

%Qy=b lösen. Q ist orthogonal Q^(-1)=Q'
y=Q'*b;

%Rx=y durch Rückwärtssubstitution lösen.
x(n,1)=y(n,1)/R(n,n);
sum=0;
for i=(n-1):-1:1
for j=i+1:n
sum=sum+R(i,j)*x(j,1);
end
x(i,1)=(y(i,1)-sum)/R(i,i);
sum=0;
end
x

%Kontrolle
xm=Ab
Allerdings muss er noch Fehler enthalten, denn ich sehe zum Einen unseren Lösungsvektor noch nicht, und zum anderen funktioniert der Ablauf nach der Frage, ob noch ein zweites LGS gelöst werden soll, nicht richtig: Gibt man dort nämlich 1 (also ja) ein, so passiert nichts (genauso bei 0)..woran liegt das?Beste Grüsse, tausend Dank und eine gute Nacht!
tigerbine Auf diesen Beitrag antworten »
RE: Lineares Gleichungssystem
Bei mir läuft es. Bei dir kommt das "end" zu früh.

kompletter coder

code:
1:
2:
3:
4:
5:
6:
7:
8:
9:
10:
11:
12:
13:
14:
15:
16:
17:
18:
19:
20:
21:
22:
23:
24:
25:
26:
27:
28:
29:
30:
31:
32:
33:
34:
35:
36:
37:
38:
39:
40:
41:
42:
43:
44:
45:
46:
47:
48:
49:
50:
51:
52:
53:
54:
55:
56:
57:
58:
59:
60:
61:
62:
63:
64:
65:
66:
67:
68:
69:
70:
71:
72:
73:
74:
75:
76:
77:
78:
79:
80:
81:
82:
83:
84:
85:
86:
87:
88:
89:
90:
91:
92:
93:
94:
95:
96:
97:
98:
%QR-Zerlegung mit Householder

clear all;
disp('Es wird eine QR-Zerlegung mit Householder berechnet.')
disp('====================================================')
disp('---------------------------')
disp('| H=I-ß*uu^T              |')
disp('| u=x+sign(x1)*||x||_2*e1 |')
disp('| ß=2/(u^Tu)              |')
disp('---------------------------')
disp(' ')
A0=input('Matrix A eingeben: ')
disp(' ')
n=length(A0);
m=length(A0');
A=A0;
min=min(m-1,n);
Q=eye(n);

for k=1:min
    fprintf('Durchgang %g \n', k)
       disp('================')
       
       clear x;
       clear u;
       clear beta;
       clear e1;
       
       e1=zeros(m+1-k,1);
       for i=k:min+1
           x(i+1-k,1)=A(i,k);
       end
       x
       x2=norm(x,2);
       signx=sign(x(1,1))
       e1(1,1)=1;
       e1
       
       u=x+signx.*x2.*e1
       beta=2/(u'*u)
       
       id=eye(n+1-k);
       disp(' ')
       H=id-beta.*(u*u')
              
       Qk=eye(n);
       for i=k:n
           for j=k:n
               Qk(i,j)=H(i+1-k,j+1-k);
           end
       end
       Qk
       Q=Qk*Q;
       Ak=Qk*A
       A=Ak;
       disp(' ')
       pause       
 end
 R=A;
 A=A0;
 Q=Q';
 disp(' ')
 disp('Die QR-Zerlegung mit Householder')
 disp('--------------------------------')
 A
 Q
 R
 
 %Berechnung einer LGS Lösung
 disp(' ');
 ans=input('Soll noch ein LGS gelöst werden? 0-nein, 1-ja: ');
 
 if ans==0
     break
 else
 b=input('Vektor b eingeben: ')
 
 %Qy=b lösen. Q ist orthogonal Q^(-1)=Q'
   y=Q'*b;
 
 %Rx=y durch Rückwärtssubstitution lösen.
 x(n,1)=y(n,1)/R(n,n);
 sum=0;
 for i=(n-1):-1:1
     for j=i+1:n         
     sum=sum+R(i,j)*x(j,1);
     end
     x(i,1)=(y(i,1)-sum)/R(i,i);
     sum=0;
 end
 x
 
 %Kontrolle
 xm=A\b
        
 end
 
Manuel20 Auf diesen Beitrag antworten »
RE: Lineares Gleichungssystem
Hmm...das kann aber nicht dein ganzer Code sein, oder?
..wenn ich den verwende, kann ich am Anfang zwar eine Matrix A eingeben, aber der Vektor b nicht.
Mein Code sieht nun folgendermassen aus:
code:
1:
2:
3:
4:
5:
6:
7:
8:
9:
10:
11:
12:
13:
14:
15:
16:
17:
18:
19:
20:
21:
22:
23:
24:
25:
26:
27:
28:
29:
30:
31:
32:
33:
34:
35:
36:
37:
38:
39:
40:
41:
42:
43:
44:
45:
46:
47:
48:
49:
50:
51:
52:
53:
54:
55:
56:
57:
58:
59:
60:
61:
62:
63:
64:
65:
66:
67:
68:
69:
70:
71:
72:
73:
74:
75:
76:
77:
78:
79:
80:
81:
82:
83:
84:
85:
86:
87:
88:
89:
90:
91:
92:
93:
94:
95:
96:
97:
98:
99:
100:
101:
102:
103:
104:
105:
106:
107:
108:
clear all;
disp('Es wird eine QR-Zerlegung mit Householder berechnet.')
disp('====================================================')
disp('---------------------------')
disp('| H=I-ß*uu^T |')
disp('| u=x+sign(x1)*||x||_2*e1 |')
disp('| ß=2/(u^Tu) |')
disp('---------------------------')
disp(' ')
A0=input('Matrix A eingeben: ')
b=input('Vektor b eingeben: ');
disp(' ')
disp(' ')
n=length(A0);
m=length(A0');
A=A0;
min=min(m-1,n);
Q=eye(n);

% Ist b ein Spaltenvektor ?
[q,p] = size(b);
if q~=n
error('b hat die falsche Dimension!');
end


for k=1:min
fprintf('Durchgang %g n', k)
disp('================')

clear x;
clear u;
clear beta;
clear e1;

e1=zeros(m+1-k,1);
for i=k:min+1
x(i+1-k,1)=A(i,k);
end
x
x2=norm(x,2);
signx=sign(x(1,1))
e1(1,1)=1;
e1

u=x+signx.*x2.*e1
beta=2/(u'*u)

id=eye(n+1-k);
disp(' ')
H=id-beta.*(u*u')

Qk=eye(n);
for i=k:n
for j=k:n
Qk(i,j)=H(i+1-k,j+1-k);
end
end
Qk
Q=Qk*Q;
Ak=Qk*A
A=Ak;
disp(' ')
pause
end


R=A;
A=A0;
Q=Q';
disp(' ')
disp('Die QR-Zerlegung mit Householder')
disp('--------------------------------')
A
Q
R


%Berechnung einer LGS Lösung
disp(' ');
ans=input('Soll noch ein LGS gelöst werden? 0-nein, 1-ja: ');

if ans==0
break
else
b=input('Vektor b eingeben: ')

%Qy=b lösen. Q ist orthogonal Q^(-1)=Q'
y=Q'*b;

%Rx=y durch Rückwärtssubstitution lösen.
x(n,1)=y(n,1)/R(n,n);
sum=0;
for i=(n-1):-1:1
for j=i+1:n
sum=sum+R(i,j)*x(j,1);
end
x(i,1)=(y(i,1)-sum)/R(i,i);
sum=0;
end
x

%Kontrolle
xm=Ab

end

Es stimmt aber nach wie vor etwas nicht, denn irgendwie erreicht der Code den ganzen LGS-Lösungs-Abschnitt nicht..
Woran liegt das? Die end's sollten nun eigentlich stimmen..
tigerbine Auf diesen Beitrag antworten »

Mein Programm lauft durch bei mir.

code:
1:
2:
3:
4:
5:
6:
7:
8:
9:
10:
11:
12:
13:
14:
15:
16:
17:
18:
19:
20:
21:
22:
23:
24:
25:
26:
27:
28:
29:
30:
31:
32:
33:
34:
35:
36:
37:
38:
39:
40:
41:
42:
43:
44:
45:
46:
47:
48:
49:
50:
51:
52:
53:
54:
55:
56:
57:
58:
59:
60:
61:
62:
63:
64:
65:
66:
67:
68:
69:
70:
71:
72:
73:
74:
75:
76:
77:
78:
79:
80:
81:
82:
83:
84:
85:
86:
87:
88:
89:
90:
91:
92:
93:
94:
95:
96:
97:
98:
99:
100:
101:
102:
103:
104:
105:
106:
107:
108:
109:
110:
111:
112:
113:
114:
115:
116:
117:
118:
119:
120:
121:
122:
123:
124:
125:
126:
127:
128:
129:
130:
131:
132:
133:
134:
135:
136:
137:
138:
139:
140:
141:
142:
143:
144:
145:
146:
147:
148:
149:
150:
151:
152:
153:
154:
155:
156:
157:
158:
159:
160:
161:
162:
163:
164:
165:
166:
167:
168:
169:
170:
171:
172:
173:
174:
 QRmitH
Es wird eine QR-Zerlegung mit Householder berechnet.
====================================================
---------------------------
| H=I-ß*uu^T              |
| u=x+sign(x1)*||x||_2*e1 |
| ß=2/(u^Tu)              |
---------------------------
 
Matrix A eingeben: [-3,2,0;4,-1,5;0,1,-1]

A0 =

    -3     2     0
     4    -1     5
     0     1    -1

 
Durchgang 1 
================

x =

    -3
     4
     0


signx =

    -1


e1 =

     1
     0
     0


u =

    -8
     4
     0


beta =

   0.02500000000000

 

H =

  -0.60000000000000   0.80000000000000                  0
   0.80000000000000   0.60000000000000                  0
                  0                  0   1.00000000000000


Qk =

  -0.60000000000000   0.80000000000000                  0
   0.80000000000000   0.60000000000000                  0
                  0                  0   1.00000000000000


Ak =

   5.00000000000000  -2.00000000000000   4.00000000000000
  -0.00000000000000   1.00000000000000   3.00000000000000
                  0   1.00000000000000  -1.00000000000000

 
Durchgang 2 
================

x =

     1
     1


signx =

     1


e1 =

     1
     0


u =

   2.41421356237309
   1.00000000000000


beta =

   0.29289321881345

 

H =

  -0.70710678118655  -0.70710678118655
  -0.70710678118655   0.70710678118655


Qk =

   1.00000000000000                  0                  0
                  0  -0.70710678118655  -0.70710678118655
                  0  -0.70710678118655   0.70710678118655


Ak =

   5.00000000000000  -2.00000000000000   4.00000000000000
   0.00000000000000  -1.41421356237309  -1.41421356237309
   0.00000000000000                  0  -2.82842712474619

 
 
Die QR-Zerlegung mit Householder
--------------------------------

A =

    -3     2     0
     4    -1     5
     0     1    -1


Q =

  -0.60000000000000  -0.56568542494924  -0.56568542494924
   0.80000000000000  -0.42426406871193  -0.42426406871193
                  0  -0.70710678118655   0.70710678118655


R =

   5.00000000000000  -2.00000000000000   4.00000000000000
   0.00000000000000  -1.41421356237309  -1.41421356237309
   0.00000000000000                  0  -2.82842712474619

 
Soll noch ein LGS gelöst werden? 0-nein, 1-ja: 1
Vektor b eingeben: [1;2;3]

b =

     1
     2
     3


x =

   1.50000000000000
   2.75000000000000
  -0.25000000000000


xm =

   1.50000000000000
   2.75000000000000
  -0.25000000000000


Bei dir steht im Test xm=Ab, da muss doch xm=A\b stehen. Das ist doch nur der Test, was mit dem Löser "\" raus kommen würde. Es läuft aber bis dahin durch.
Manuel20 Auf diesen Beitrag antworten »

Super - es hat tatsächlich an dem einen "\" gelegen smile
Jupiiiiiiiii Tanzen

Eine Frage noch zum Schluss:
Wenn ich das Programm am folgenden Beispiel testen möchte:



Wie gebe ich das am besten ein?
Der "Sinn" liegt darin, dass man dann die maximale Dimension n herausfinden soll, für welche die QR Zerlegung die Lösung von Ax=b in weniger als einer Minute liefert.. smile
tigerbine Auf diesen Beitrag antworten »

Blockiere die manuelle Eingabe mit %

blockiere die Textausgaben mit; zeige nur das Ergebnis an (ohne A\b berechnen zu lassen).

schreibe drum rum eine Schleife, die die MAtrix und den Vektor passend erzeugt.

Wie man nun Zeit als abbruch nimmt, weiß ich nciht. Würde eben beginn und ende der Schleife anpassen.
Manuel20 Auf diesen Beitrag antworten »

Ist es nicht auch möglich, das Programm so zu verändern, dass es im Hintergrund automatisch arbeitet (sprich: die Dimension selber immer erhöht), und dass man dann manuell die Zeit (hier eben nach 1 min) abbricht, zB mit der Eingabe einer 0 oder was auch immer...?
tigerbine Auf diesen Beitrag antworten »

ich sagte doch, schleife. mit strg+c kannst du abbrechen.
Manuel20 Auf diesen Beitrag antworten »

Oke, ich habe allerdings noch 2, 3 Fragen:

Wie definiere ich die Matrix A richtig?
Ich habe es mal so gemacht, aber das "e" existiert so in matlab nicht.
(auch mit "eul" funktioniert das Ganze nicht wirklich)
A0=[e^(-abs(i-j))];

Dann zum Vektor:
Versteht matlab folgendes?:
b=[1]^T
oder wie müsste ich das implementieren?
(gemeint ist das: b=(1,1,......,1)^T )

Wieso kann ich nicht einfach vor der Definition von A folgendes schreiben?:
for i<=1
for j<=n
(...)
end
end

Besten Dank für die Hilfe!
tigerbine Auf diesen Beitrag antworten »

exp heißt die Funktion. Bastel damit was. Bitte bemühe auch mal das Internet für reine Fragen zur matlab Sprache. Mit ones kann man einen Vektr/Matrix nur aus 1er erzeugen.
Manuel20 Auf diesen Beitrag antworten »

Hey!

Danke für all Deine Tipps und Erklärungen!

Könntest du Dir mein Code (hoffentlich ein letztes Mal) nochmals anschauen?
Irgendwas ist (evtl gar grundlegend) noch falsch, da bei mir nach 2sec die Fehlermeldung:
??? Attempt to execute SCRIPT ans as a function: erscheint.

code:
1:
2:
3:
4:
5:
6:
7:
8:
9:
10:
11:
12:
13:
14:
15:
16:
17:
18:
19:
20:
21:
22:
23:
24:
25:
26:
27:
28:
29:
30:
31:
32:
33:
34:
35:
36:
37:
38:
39:
40:
41:
42:
43:
44:
45:
46:
47:
48:
49:
50:
51:
52:
53:
54:
55:
56:
57:
58:
59:
60:
61:
62:
63:
64:
65:
66:
67:
68:
69:
70:
71:
72:
73:
74:
75:
76:
77:
78:
79:
80:
81:
82:
83:
84:
85:
86:
87:
88:
89:
90:
91:
92:
93:
94:
95:
96:
97:
98:
99:
100:
101:
102:
103:
104:
105:
106:
107:
108:
109:
110:
111:
112:
113:
114:
115:
116:
%QR-Zerlegung mit Householder

clear all;
disp('Es wird eine QR-Zerlegung mit Householder berechnet.')
disp('====================================================')
disp('---------------------------')
disp('| H=I-ß*uu^T |')
disp('| u=x+sign(x1)*||x||_2*e1 |')
disp('| ß=2/(u^Tu) |')
disp('---------------------------')
disp(' ')

n=inf;
for i=1:n
for j=1:n
A0=[exp(-abs(i-j))];

b=[ones];
disp(' ')
disp(' ')
n=length(A0);
m=length(A0');
A=A0;
min=min(m-1,n);
Q=eye(n);

% Ist b ein Spaltenvektor ?
[q,p] = size(b);
if q~=n
error('b hat die falsche Dimension!');
end


for k=1:min
fprintf('Durchgang %g n', k)
disp('================')

clear x;
clear u;
clear beta;
clear e1;

e1=zeros(m+1-k,1);
for i=k:min+1
x(i+1-k,1)=A(i,k);
end
x
x2=norm(x,2);
signx=sign(x(1,1))
e1(1,1)=1;
e1

u=x+signx.*x2.*e1
beta=2/(u'*u)

id=eye(n+1-k);
disp(' ')
H=id-beta.*(u*u')

Qk=eye(n);
for i=k:n
for j=k:n
Qk(i,j)=H(i+1-k,j+1-k);
end
end
Qk
Q=Qk*Q;
Ak=Qk*A
A=Ak;
disp(' ')
pause
end


R=A;
A=A0;
Q=Q';
disp(' ')
disp('Die QR-Zerlegung mit Householder')
disp('--------------------------------')
A
Q
R


%Berechnung einer LGS Lösung
disp(' ');
%ans=input('Soll noch ein LGS gelöst werden? 0-nein, 1-ja: ');

if ans==0
break
else
%b=input('Vektor b eingeben: ')

%Qy=b lösen. Q ist orthogonal Q^(-1)=Q'
y=Q'*b;

%Rx=y durch Rückwärtssubstitution lösen.
x(n,1)=y(n,1)/R(n,n);
sum=0;
for i=(n-1):-1:1
for j=i+1:n
sum=sum+R(i,j)*x(j,1);
end
x(i,1)=(y(i,1)-sum)/R(i,i);
sum=0;
end
x

%Kontrolle
%xm=Ab

end
end
end
Neue Frage »
Antworten »



Verwandte Themen

Die Beliebtesten »
Die Größten »
Die Neuesten »