| 18.03.2010, 23:15 |
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! |
| 19.03.2010, 09:56 |
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. 
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
|
|
|
| 19.03.2010, 12:39 |
Manuel20 |
Auf diesen Beitrag antworten » |
RE: Lineares Gleichungssystem
Hallo tigerbine!
Ja, dass die QR-Zerlegung wunderbar klappt, habe ich gesehen - super! 
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? |
| 19.03.2010, 12:54 |
tigerbine |
Auf diesen Beitrag antworten » |
RE: Lineares Gleichungssystem
Wann ist egal. Solange es geschieht, bevor du drauf zugreifst. 
| code: |
1:
2:
3:
|
b=input('Vektor b eingeben: ');
|
|
Dann musst du eben nur überlegen, wie du aus QRx=b x berechnen kannst. |
| 19.03.2010, 14:12 |
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
..wuerde das klappen? (bzw wie (codemaessig)) |
| 19.03.2010, 14:21 |
tigerbine |
Auf diesen Beitrag antworten » |
RE: Lineares Gleichungssystem
Durch eine Matrix teilen.... aua.... 
Die erste Idee ist die richtige. Und ich sagte das ja auch schon.
| Zitat: |
| Q ist ja wunderschön.
Und das letzte LGS dann, wegen der Gestalt von R durch Rückwärtssubstitution. |
|
| Anzeige | |
|
|
| 19.03.2010, 19:12 |
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.. |
| 19.03.2010, 19:24 |
tigerbine |
Auf diesen Beitrag antworten » |
RE: Lineares Gleichungssystem
Verstehe das Problem nicht.
Was soll das sein
Ich sehe nirgends das doppelte Lösen eines LGS. Dann hättest du ja auch x=A\b eingeben können. |
| 19.03.2010, 21:06 |
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? |
| 19.03.2010, 21:18 |
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.
Das mit A\b war als Witz gemeint. Dann braucht man auch nicht doe QR-Zerlegung zu berechnen, wenn man den Matlablöser "\" nimmt.
|
| 19.03.2010, 22:32 |
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! |
| 19.03.2010, 22:35 |
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.
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. |
| 19.03.2010, 22:45 |
Manuel20 |
Auf diesen Beitrag antworten » |
RE: Lineares Gleichungssystem
Hehe
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. |
| 19.03.2010, 22:51 |
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. |
| 19.03.2010, 22:59 |
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] |
| 19.03.2010, 23:03 |
tigerbine |
Auf diesen Beitrag antworten » |
RE: Lineares Gleichungssystem
Gib mal nur die Matrix ein und dann rank(A). Schlechtes Beispiel, was du genommen hast.
|
| 19.03.2010, 23:19 |
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... |
| 19.03.2010, 23:22 |
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 |
| 19.03.2010, 23:30 |
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. |
| 19.03.2010, 23:43 |
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. |
| 20.03.2010, 18:16 |
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? |
| 20.03.2010, 19:04 |
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
stehen hast, gibt man x eh schon aus. Du kannst das optisch so aufblähen, wie du magst.
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. 
|
| 20.03.2010, 21:28 |
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
|
| 20.03.2010, 22:13 |
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. |
| 20.03.2010, 23:10 |
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.. 
|
| 20.03.2010, 23:29 |
tigerbine |
Auf diesen Beitrag antworten » |
RE: Lineares Gleichungssystem
Ich versteh einfach nicht, was du nicht verstehst. 
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'. |
| 21.03.2010, 00:39 |
Manuel20 |
Auf diesen Beitrag antworten » |
RE: Lineares Gleichungssystem
ahhhh..jetzt ist vieles klarer, zumindest was das Verständnis angeht
Weil Q orthogonal ist, lassen sich die Inverse einfach berechnen, was dann die (komplex) transponierte Matrix der Matrix Q ist
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. |
| 21.03.2010, 00:42 |
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(' ')
|
|
|
| 21.03.2010, 00:44 |
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. |
| 21.03.2010, 01:32 |
Manuel20 |
Auf diesen Beitrag antworten » |
RE: Lineares Gleichungssystem
Wow, du bist ein echtes Genie!
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! |
| 21.03.2010, 01:41 |
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
|
|
|
| 21.03.2010, 11:08 |
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.. |
| 21.03.2010, 12:08 |
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. |
| 21.03.2010, 17:28 |
Manuel20 |
Auf diesen Beitrag antworten » |
Super - es hat tatsächlich an dem einen "\" gelegen
Jupiiiiiiiii 
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..
|
| 21.03.2010, 17:32 |
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. |
| 21.03.2010, 18:07 |
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...? |
| 21.03.2010, 18:08 |
tigerbine |
Auf diesen Beitrag antworten » |
ich sagte doch, schleife. mit strg+c kannst du abbrechen. |
| 21.03.2010, 19:46 |
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! |
| 21.03.2010, 20:43 |
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. |
| 21.03.2010, 22:16 |
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 |
|
|
