Octave

      So that's going to be my own thread where I post stuff about,and my learning status of, Octave.
      No one's going to read this anyway so who cares!
      What is Octave? Octave is an interpreting language!

      Octave get it here

      It's easy to learn and simple to use!

      Basic Tutorial

      What is the purpose of this thread may the lone reader think: Well I don't know it's a good way to train my english and to keep track of my progress.
      Also Pyro might read this and give me some help lululolol








      Come along, follow me, as I lead through the darkness.

      Post was edited 1 time, last by “Sparrow” ().

      Since I have a broken leg and lots of time,
      Here is the first update:

      I wrote my first 2 functions.
      Back and forward substitution.

      forward substitution:

      Inputs: A Lower triangular matrix and a column vector

      Source Code

      1. function x = vorsub(R,z)
      2. [a,b] = size(R);
      3. [c,d] = size(z);
      4. if a != b,
      5. error("Es wird eine n x n Matrix erwartet" )
      6. end;
      7. if a != c,
      8. error("R und z haben nicht die gleiche Dimension")
      9. end;
      10. if d != 1,
      11. error("Es wird ein n x 1 Vektor erwartet" )
      12. end;
      13. x = zeros(a,1);
      14. x(1) = z(1) / R(1,1);
      15. for i=2:a,
      16. x(i) =1/R(i,i) * (z(i) - sum(R(i,i:a-1) * x(i-1)));
      17. end;
      18. endfunction


      backward substitution:

      Inputs: A Triangular matrix and a column vector

      Source Code

      1. function x = ruecksub(R,z)
      2. [a,b] = size(R);
      3. [c,d] = size(z);
      4. if a != b,
      5. error("Es wird eine n x n Matrix erwartet" )
      6. end;
      7. if a != c,
      8. error("R und z haben nicht die gleiche Dimension")
      9. end;
      10. if d != 1,
      11. error("Es wird ein n x 1 Vektor erwartet" )
      12. end;
      13. x = zeros(a,1);
      14. x(a) = z(a)/R(a,a);
      15. for i = a-1:-1:1,
      16. x(i) = 1/R(i,i) * (z(i) - sum(R(i,i+1:a) * x(i+1)));
      17. end;
      18. endfunction
      Come along, follow me, as I lead through the darkness.

      Post was edited 1 time, last by “Sparrow” ().

      Update:

      LU decomposition without pivoting:

      Input:
      n x n matrix

      Source Code

      1. function [L,A] = LRzerl (A)
      2. % input ist eine n x n Matrix, errechnet wird die Linke Seite in einen
      3. % neuen Vector während die Matrix A mit R ueberschrieben wird.
      4. % => Ausgabe ist L = Linksseitige Matrix, A = Rechtsseitige Matrix
      5. [a,b] = size(A);
      6. %ueberpruefung ob n x n matrix
      7. if a != b,
      8. error ("Es wird eine nxn Matrix erwartet");
      9. end;
      10. % Identitätsmatrix n x n
      11. L = eye(a);
      12. for k = 1:a-1,
      13. for i = k+1:a,
      14. L(i,k) = A(i,k)/A(k,k);
      15. for j = i-1:a,
      16. % R ueberschreibt A
      17. A(i,j) = A(i,j)-L(i,k) * A(k,j);
      18. end;
      19. end;
      20. end;
      21. endfunction
      Come along, follow me, as I lead through the darkness.
      Forget what I wrote up there, this is the new stuff and it's better!


      Source Code

      1. function [L,R] = LRzerl (A)
      2. % input ist eine n x n Matrix, errechnet wird die Linke Seite in einen
      3. % neuen Vector während die Matrix A mit R ueberschrieben wird.
      4. % => Ausgabe ist L = Linksseitige Matrix, A = Rechtsseitige Matrix
      5. [a,b] = size(A);
      6. %ueberpruefung ob n x n matrix
      7. if a != b,
      8. error ("Es wird eine nxn Matrix erwartet");
      9. end;
      10. % Identitätsmatrix n x n
      11. L = eye(a);
      12. %LR Zerlegung nach Pseudocode von Num1Skript.pdf Seite 14
      13. for k = 1:a-1,
      14. for i = k+1:a,
      15. L(i,k) = A(i,k)/A(k,k);
      16. for j = k+1:a,
      17. % R ueberschreibt A, hier entstehen die Werte R in A, der Rest der Werte
      18. % in A links unten bleiben gleich
      19. A(i,j) = A(i,j)-L(i,k) * A(k,j);
      20. end;
      21. end;
      22. end;
      23. %links unten null setzen
      24. R = triu(A);
      25. endfunction
      Come along, follow me, as I lead through the darkness.
      and with pivoting...:


      Source Code

      1. function [L,R,P] = pivoLRzerl (A)
      2. %input ist eine n x n Matrix, errechnet wird die Linke Seite in einen
      3. % neuen Vector während R nur A kopieren und unten links null setzen muss
      4. % der Pivotvektor hat anfangs die Zahlen 1:a und 0wird innerhalb der funktion
      5. % bestimmt und später überall dort, wo berechnet mit 1 besetzt
      6. % => Ausgabe ist L = Linksseitige Matrix, R = Rechtsseitige Matrix, P = Pivotvektor
      7. [a,b] = size(A);
      8. %ueberpruefung ob n x n matrix
      9. if a != b,
      10. error ("Es wird eine nxn Matrix erwartet");
      11. end;
      12. %Vektor 1xn erstellen
      13. p = 1:a;
      14. %Identitätsmatrix n x n
      15. L = eye(a);
      16. %LR Zerlegung nach Pseudocode von Num1Skript.pdf Seite 14
      17. for k = 1:a-1,
      18. %Vertauschung nach "Matlab-Notation" in Num1Skript.pdf Seite 19
      19. %bestimme den absolut groessten Wert in Zeile k und tausche ihn dann
      20. [c,d] = max(abs(A(p(k:end),k)));
      21. %vertauschen
      22. ps = p(k);
      23. p(k) = p(k-1+d);
      24. p(k-1+d) = ps;
      25. %mit vertauschten Spalten LR Zerlegung fortfuehren
      26. for i = k+1:a,
      27. A(p(i),k) = A(p(i),k)/A(p(k),k);
      28. for j = k+1:a,
      29. A(p(i),j) = A(p(i),j)-A(p(i),k) * A(p(k),j);
      30. end;
      31. end;
      32. end;
      33. %L kopieren und mittig einsen einfuegen
      34. L = L + tril(A(p,1:a),-1);
      35. %links unten null setzen
      36. R = triu(A(p,1:a));
      37. %P berechnen
      38. P = eye(a); P(1:a,p)=P;
      39. endfunction
      Come along, follow me, as I lead through the darkness.
      Ok so I decided to completly re-organize my project!
      Why?
      Because it all comes down to the Gaussian elimination with or without pivoting and in that case it's more comfortable to have it all running with each other and having the pivoting stuff as an secondary option inside of the function.
      Also the output is simplified by just giving you the result and the matrices are now saved compact inside of A.

      forward substitution:

      Source Code

      1. function z = vorsub(L,b)
      2. % Input ist R = n x n Matrix, z = n x 1 Vektor
      3. % Output x Vektor mit n x 1 Loesungen
      4. [a,f] = size(L);
      5. [c,d]= size(b);
      6. if a != f,
      7. error("Es wird eine n x n Matrix erwartet" )
      8. end;
      9. if a != c,
      10. error("R und z haben nicht die gleiche Dimension")
      11. end;
      12. if d != 1,
      13. error("Es wird ein n x 1 Vektor erwartet" )
      14. end;
      15. %nullvektor mit der richtigen groese
      16. z = zeros(a,1);
      17. z(1) = b(1);
      18. %nach dem Algorithmus der Vorwaertssubstition
      19. for i=2:1:a,
      20. z(i) =(b(i) - L(i,1:i-1) * z(1:i-1));
      21. end;
      22. endfunction


      backward substitution:

      Source Code

      1. function z = ruecksub(R,x)
      2. % Input ist R = n x n Matrix, z = n x 1 Vektor
      3. % Output x Vektor mit n x 1 Loesungen
      4. [a,b] = size(R);
      5. [c,d] = size(x);
      6. if a != b,
      7. error("Es wird eine n x n Matrix erwartet" )
      8. end;
      9. if a != c,
      10. error("R und z haben nicht die gleiche Dimension")
      11. end;
      12. if d != 1,
      13. error("Es wird ein n x 1 Vektor erwartet" )
      14. end;
      15. %nullvektor mit der richtigen groese
      16. z = zeros(a,1);
      17. if abs(R(a,a)) < eps,
      18. error("Nicht durchfuehrbar");
      19. end;
      20. z(a) = x(a)/R(a,a);
      21. %nach dem Algorithmus der Rueckwaertssubstition
      22. for i = a-1:-1:1,
      23. if abs(R(i,i)) < eps,
      24. error("Nicht durchfuehrbar");
      25. end;
      26. z(i) = 1/R(i,i) * (x(i) - R(i,i+1:a) * z(i+1:a));
      27. end;
      28. endfunction


      LU decomposition with or without pivoting

      Source Code

      1. function [A,p] = pivoLRzerl(A,s)
      2. % input ist eine n x n Matrix, errechnet wird die Linke Seite in einen
      3. % neuen Vector während R nur A kopieren und unten links null setzen muss
      4. % der Pivotvektor hat anfangs die Zahlen 1:a und 0wird innerhalb der funktion
      5. % bestimmt und später überall dort, wo berechnet mit 1 besetzt
      6. % => Ausgabe ist L = Linksseitige Matrix, R = Rechtsseitige Matrix, P = Pivotvektor
      7. [a,b] = size(A);
      8. %ueberpruefung ob n x n matrix
      9. if a != b,
      10. error ("Es wird eine nxn Matrix erwartet");
      11. end;
      12. %Vektor 1xn erstellen
      13. p = (1:a)';
      14. %Identitätsmatrix n x n
      15. %LR Zerlegung nach Pseudocode von Num1Skript.pdf Seite 14
      16. for k = 1:a-1,
      17. %Pivot Ja - Nein?
      18. % 1 == Ja
      19. if s == 1,
      20. %Vertauschung nach "Matlab-Notation" in Num1Skript.pdf Seite 19
      21. %bestimme den absolut groessten Wert in Zeile k und tausche ihn dann
      22. [c,d] = max(abs(A(p(k:end),k)));
      23. %vertauschen
      24. ps = p(k);
      25. p(k) = p(k-1+d);
      26. p(k-1+d) = ps;
      27. end;
      28. %mit vertauschten Spalten LR Zerlegung fortfuehren
      29. for i = k+1:a,
      30. A(p(i),k) = A(p(i),k)/A(p(k),k);
      31. for j = k+1:a,
      32. A(p(i),j) = A(p(i),j)-A(p(i),k) * A(p(k),j);
      33. end;
      34. end;
      35. end;
      36. endfunction


      Gaussian Elimination with or without pivoting

      Source Code

      1. function x = gauss(A,b,s)
      2. a = size(A,1);
      3. %LR zerlegung mit Pivo?...
      4. if (s == 1),
      5. %zerlegen mit Pivot
      6. [A,p] = pivoLRzerl(A,1);
      7. e = b(p);
      8. %vorwärtssubstitution
      9. z = vorsub(A,e);
      10. else,
      11. %...ansonsten ohne
      12. %gleich nur ohne Pivo (daher 2)
      13. [A,p] = pivoLRzerl(A,2);
      14. z = vorsub(A,b);
      15. end;
      16. %rueckwaertssubstitution
      17. x = ruecksub(A,z);
      18. endfunction
      Come along, follow me, as I lead through the darkness.
      Update:

      The Gauss-Seidel and SOR method.

      Inputs: Matrix a, right side b, start value x0, tolerance tol, maximum iterations maxiter, Over-relaxation parameter w
      Output: Solution vector x, iteration amount g

      Source Code

      1. function [x,g] = seidel (A,b,x0,tol,maxiter,w)
      2. [a,s] = size(A);
      3. x = x0;
      4. r = norm(b - A*x);
      5. g = 0;
      6. while (r > tol && g < maxiter ),
      7. g = g + 1;
      8. for j = 1:a,
      9. x(j) = (1-w)*x(j) + w/A(j,j) * (b(j) - A(j,1:j-1) * x(1:j-1) - A(j,j+1:a) * x(j+1:a));
      10. end;
      11. r = norm(b-A*x);
      12. end;
      Come along, follow me, as I lead through the darkness.
      The Poisson equation plotted.
      Minimum iterations search solved with the SOR-Function and over the spectral radius, also plotted.

      Source Code

      1. N=10;
      2. h=1/(N+1);
      3. x=0:h:1;
      4. y=0:h:1;
      5. [X,Y] = meshgrid(x,y);
      6. TN=sparse(4*eye(N));
      7. for k=2:N
      8. TN(k-1,k)=-1;
      9. TN(k,k-1)=-1;
      10. end
      11. IN=sparse(eye(N));
      12. A=[TN];
      13. for k=2:N
      14. A(N*(k-2)+1:N*(k-1),N*(k-1)+1:N*k)=-IN;
      15. A(N*(k-1)+1:N*k,N*(k-2)+1:N*k )= [-IN TN];
      16. end
      17. % Rechte Seite für die PDE u_xx + u_yy = 5
      18. f=-h^2*5*ones(N^2,1);
      19. % Einbeziehung der Randwerte:
      20. % u(x, y)=0;
      21. u0i=zeros(N+2,1);
      22. uN1i=zeros(N+2,1);
      23. ui0=zeros(N+2,1);
      24. uiN1=zeros(N+2,1);
      25. b=zeros(N^2,1);
      26. b(1:N)=ui0(2:N+1);
      27. b(1)=b(1)+u0i(2);
      28. b(N)=b(N)+uN1i(2);
      29. for k=2:N-1
      30. b((k-1)*N+1)=u0i(k+1);
      31. b(k*N)=uN1i(k+1);
      32. end
      33. b(N^2-N+1:N^2)=uiN1(2:N+1);
      34. b(N^2-N+1)=b(N^2-N+1)+u0i(N+1);
      35. b(N^2)=b(N^2)+uN1i(N+1);
      36. c=b+f;
      37. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
      38. % mit MATLAB-Solver
      39. %sol=A\c;
      40. %Aufgabe 4
      41. sol = seidel(A,c,b,1.0*e^-3,50,1.299);
      42. %Aufgabe 5
      43. axis([0 2 0 25]);
      44. hold on;
      45. for w0 = 0.1:0.1:1.9,
      46. [z,y] = seidel(A,c,b,1.0*e^-3,20,w0);
      47. plot(w0,y,"*","markersize",8);
      48. hold on;
      49. end;
      50. title('Aufgabe 5')
      51. xlabel('W')
      52. ylabel('maxiter')
      53. grid on
      54. figure
      55. %Aufgabe 6
      56. %Nach Aufgabenstellung
      57. D=diag(diag(A));
      58. L=tril(A)-D;
      59. R=triu(A)-D;
      60. [l,k] = size(A);
      61. I = eye(l,k);
      62. axis([0 2 0 1.5]);
      63. hold on;
      64. for w0 = 0.1:0.1:1.9,
      65. B = inv(I + w0*inv(D)*L) * ((1-w0)*I-w0*inv(D)*R);
      66. spek = max(abs((eig(B))));
      67. plot(w0,spek,"*","markersize",8);
      68. hold on;
      69. end;
      70. title('Aufgabe 6')
      71. xlabel('W')
      72. ylabel('spek')
      73. grid on
      74. figure
      75. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
      76. Z=sol(1:N)';
      77. for k=2:N
      78. Z=[Z ; sol((k-1)*N+1:k*N)'];
      79. end
      80. Z=[ui0';
      81. u0i(2:N+1),Z, uN1i(2:N+1);
      82. uiN1'];
      83. surfc(X,Y,Z)
      84. colormap(jet)
      85. xlabel('x-Achse')
      86. ylabel('y-Achse')
      87. zlabel('u-Achse')
      Come along, follow me, as I lead through the darkness.
      Plotting could be solved with vectors but I'm not sure about plotting them as points and not in a straight line. Hard to read out the minimum value without seeing if there is an actual point or not...
      Come along, follow me, as I lead through the darkness.
      So... an example for the Peaucellier mechanism with the newton method?!
      I'm having problems giving functions to my functions, mainly because Octave is shit... thanks @Cl0wN -.-

      Brainfuck Source Code

      1. ​ function[] = Peaucellier
      2. %Global parameters
      3. global L1;
      4. global L2;
      5. global L3;
      6. %Global coordinates
      7. global Ax;
      8. global Ay;
      9. global Bx;
      10. global By;
      11. %Global variable
      12. global time;
      13. % Initialization of global variables
      14. L1=15;
      15. L2=20;
      16. L3=1/2*(4*(-2*L1-1/2*2^(1/2)*(L2^2)^(1/2))^2+2*L2^2)^(1/2);
      17. Ax=0;
      18. Ay=0;
      19. Bx=L1;
      20. By=0;
      21. time=0;
      22. %Initialization
      23. Cx0=2*L1;
      24. Cy0=0;
      25. Dx0=2*L1+sqrt(L2^2/2);
      26. Dy0=sqrt(L2^2/2);
      27. Ex0=2*L1+sqrt(L2^2/2);
      28. Ey0=-sqrt(L2^2/2);
      29. Fx0=2*L1+2*sqrt(L2^2/2);
      30. Fy0=0;
      31. x0=[Cx0 Cy0 Dx0 Dy0 Ex0 Ey0 Fx0 Fy0];
      32. %... Plot initial configuration
      33. PlotPeaucellier(x0)
      34. pause
      35. x=x0'; %Initialization for Newton
      36. for t=0:0.05:20
      37. hold off
      38. time=t; % Set global variable "time"
      39. % Solve nonlinear system of equations
      40. [x, NumIt] = newton(@fPeaucellier, @JPeaucellier, x,1.0e-12, 100);
      41. %... Plot Peaucellier
      42. clf
      43. PlotPeaucellier(x)
      44. pause(0.05);
      45. end
      46. hold off
      47. % -----------------------------------------------------------------------
      48. % Nested functions
      49. % -----------------------------------------------------------------------
      50. %
      51. function f = fPeaucellier(x)
      52. % Nonlinear equations describing the geometry of the Peaucelleir mechanism
      53. Cx=x(1);
      54. Cy=x(2);
      55. Dx=x(3);
      56. Dy=x(4);
      57. Ex=x(5);
      58. Ey=x(6);
      59. Fx=x(7);
      60. Fy=x(8);
      61. f = [ (Ax-Dx)^2+(Ay-Dy)^2- L3^2;
      62. (Ax-Ex)^2+(Ay-Ey)^2- L3^2;
      63. (Bx-Cx)^2+(By-Cy)^2- L1^2;
      64. (Fx-Ex)^2+(Fy-Ey)^2- L2^2;
      65. (Ex-Cx)^2+(Ey-Cy)^2- L2^2;
      66. (Cx-Dx)^2+(Cy-Dy)^2- L2^2;
      67. (Dx-Fx)^2+(Dy-Fy)^2- L2^2;
      68. Cy-L1*sin(pi/4*sin(time))];
      69. end
      70. % -----------------------------------------------------------------------
      71. function Jf = JPeaucellier(x)
      72. % Jacobian function of the nonlinear equations describing the geometry of the
      73. % Peaucelleir mechanism
      74. Cx=x(1);
      75. Cy=x(2);
      76. Dx=x(3);
      77. Dy=x(4);
      78. Ex=x(5);
      79. Ey=x(6);
      80. Fx=x(7);
      81. Fy=x(8);
      82. Jf =[ 0 0 -2*Ax+2*Dx -2*Ay+2*Dy 0 0 0 0
      83. 0 0 0 0 -2*Ax+2*Ex -2*Ay+2*Ey 0 0
      84. -2*Bx+2*Cx -2*By+2*Cy 0 0 0 0 0 0
      85. 0 0 0 0 -2*Fx+2*Ex -2*Fy+2*Ey 2*Fx-2*Ex 2*Fy-2*Ey
      86. -2*Ex+2*Cx -2*Ey+2*Cy 0 0 2*Ex-2*Cx 2*Ey-2*Cy 0 0
      87. 2*Cx-2*Dx 2*Cy-2*Dy -2*Cx+2*Dx -2*Cy+2*Dy 0 0 0 0
      88. 0 0 2*Dx-2*Fx 2*Dy-2*Fy 0 0 -2*Dx+2*Fx -2*Dy+2*Fy
      89. 0 1 0 0 0 0 0 0];
      90. end
      91. function [] = PlotPeaucellier(x)
      92. %Plot the Peaucellier mechanims
      93. Cx=x(1);
      94. Cy=x(2);
      95. Dx=x(3);
      96. Dy=x(4);
      97. Ex=x(5);
      98. Ey=x(6);
      99. Fx=x(7);
      100. Fy=x(8);
      101. plot([Ax Dx],[Ay Dy],'ro-','LineWidth',8)
      102. set(gca,'DataAspectRatio',[1 1 1])
      103. hold on
      104. plot((2*L1+2*sqrt(L2^2/2))*ones(21),-30:3:30,'-.m')
      105. xc=5/3*L1:0.1:2*L1;
      106. yc1=sqrt(L1^2-(xc-L1).^2);
      107. yc2=-sqrt(L1^2-(xc-L1).^2);
      108. plot(xc,yc1,'-.c')
      109. plot(xc,yc2,'-.c')
      110. plot([Ax Ex],[Ay Ey],'ro-','LineWidth',8)
      111. plot([Bx Cx],[By Cy],'bo-','LineWidth',8)
      112. plot([Cx Dx],[Cy Dy],'go-','LineWidth',8)
      113. plot([Cx Ex],[Cy Ey],'go-','LineWidth',8)
      114. plot([Dx Fx],[Dy Fy],'go-','LineWidth',8)
      115. plot([Ex Fx],[Ey Fy],'go-','LineWidth',8)
      116. plot(Ax,Ay,'ko','MarkerSize',9,'MarkerFaceColor','y')
      117. plot(Bx,By,'ok','MarkerSize',9,'MarkerFaceColor','y')
      118. plot(Cx,Cy,'ok','MarkerSize',8,'MarkerFaceColor','w')
      119. plot(Dx,Dy,'ok','MarkerSize',8,'MarkerFaceColor','w')
      120. plot(Ex,Ey,'ok','MarkerSize',8,'MarkerFaceColor','w')
      121. plot(Fx,Fy,'ok','MarkerSize',8,'MarkerFaceColor','w')
      122. title('Mechanismus von Peaucellier')
      123. shg
      124. end
      125. end
      Come along, follow me, as I lead through the darkness.
      Newton möh:

      Source Code

      1. ​function [x,Numiter] = newton (fun, Jfun, x0, tol, itermax)
      2. #Input: Funktionen Fun, Jfun, Startwert, Toleranz, Maximale Iterationen
      3. #Output: Lösungsvektor x, getätigte Iterationen Numiter
      4. Numiter=0;
      5. xvor = x0;
      6. fehler = 2*tol;
      7. r = 2*tol;
      8. while ( Numiter < itermax && fehler > tol && r > tol ),
      9. %Auswertung der Funktion und der Jacobi-Matrix:
      10. J = Jfun(xvor);
      11. f = fun(xvor);
      12. dx = gauss(J,-f,1);
      13. x=xvor+dx;
      14. fehler = norm (dx);
      15. r = norm(f);
      16. %Zaehler +1
      17. Numiter = Numiter + 1;
      18. xvor = x;
      19. end;
      20. if( Numiter == itermax ),
      21. fprintf('Konvergiert fuer die vorgegebene iterations - Grenze mit folgenden Werten nicht. %d \n', Numiter);
      22. end;
      23. if(fehler > tol),
      24. fprintf('Konvergiert mit der vorgegebenen Toleranzgrenze und folgendem Fehler nicht. %d \n', fehler);
      25. end;
      26. if(r > tol),
      27. fprintf('Konvergiert fuer das Residuum > als die Toleranz mit folgendem Residuum nicht. %d \n', r);
      28. end;
      29. endfunction
      Come along, follow me, as I lead through the darkness.