octave:2> A=sprand(5,5,0.2) A = Compressed Column Sparse (rows = 5, cols = 5, nnz = 5 [20%]) (2, 1) -> 0.059293 (2, 3) -> 0.98708 (4, 3) -> 0.84526 (1, 4) -> 0.20234 (1, 5) -> 0.78357 octave:3> A*rand(5,1) ans = 0.68209 0.96522 0.00000 0.79714 0.00000 octave:4> A(3,4) ans = Compressed Column Sparse (rows = 1, cols = 1, nnz = 0 [0%]) octave:5> B=rand(5) B = 0.569043 0.114655 0.034019 0.074413 0.598674 0.668069 0.530005 0.485480 0.774373 0.524977 0.491692 0.782651 0.097963 0.227708 0.718043 0.919250 0.818621 0.301697 0.341493 0.417917 0.303097 0.913648 0.705134 0.668763 0.081793 octave:6> whos Variables in the current scope: Attr Name Size Bytes Class ==== ==== ==== ===== ===== A 5x5 84 double B 5x5 200 double ans 1x1 8 double Total is 30 elements using 292 bytes octave:7> B=zeros(5,5) B = 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 octave:8> whos Variables in the current scope: Attr Name Size Bytes Class ==== ==== ==== ===== ===== A 5x5 84 double B 5x5 200 double ans 1x1 8 double Total is 30 elements using 292 bytes octave:9> B B = 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 octave:10> B=sparse(B) B = Compressed Column Sparse (rows = 5, cols = 5, nnz = 0 [0%]) octave:11> whos Variables in the current scope: Attr Name Size Bytes Class ==== ==== ==== ===== ===== A 5x5 84 double B 5x5 24 double ans 1x1 8 double Total is 5 elements using 116 bytes octave:12> B B = Compressed Column Sparse (rows = 5, cols = 5, nnz = 0 [0%]) octave:13> A=rand(5) A = 0.243873 0.920907 0.515177 0.030778 0.961502 0.662900 0.991719 0.583537 0.925679 0.140909 0.697899 0.520464 0.797578 0.720184 0.845110 0.586324 0.871377 0.838248 0.655750 0.325423 0.685445 0.724500 0.893872 0.961905 0.273673 octave:14> D=diag(A) D = 0.24387 0.99172 0.79758 0.65575 0.27367 octave:15> D=diag(diag(A)) D = Diagonal Matrix 0.24387 0 0 0 0 0 0.99172 0 0 0 0 0 0.79758 0 0 0 0 0 0.65575 0 0 0 0 0 0.27367 octave:16> spdiags(diag(A),0,5,5) ans = Compressed Column Sparse (rows = 5, cols = 5, nnz = 5 [20%]) (1, 1) -> 0.24387 (2, 2) -> 0.99172 (3, 3) -> 0.79758 (4, 4) -> 0.65575 (5, 5) -> 0.27367 octave:17> J_J=(speye(5)-inv(D)*A) % NOOOOOOOOOOOOOOOOOOOOOOOOOOOO J_J = 0.00000 -3.77617 -2.11248 -0.12621 -3.94263 -0.66844 0.00000 -0.58841 -0.93341 -0.14209 -0.87502 -0.65256 0.00000 -0.90296 -1.05960 -0.89413 -1.32883 -1.27830 0.00000 -0.49626 -2.50461 -2.64732 -3.26620 -3.51479 0.00000 octave:18> J_J=(speye(5)-diag(1./diag(D))*A) % NOOOOOOOOOO J_J = 0.00000 -3.77617 -2.11248 -0.12621 -3.94263 -0.66844 0.00000 -0.58841 -0.93341 -0.14209 -0.87502 -0.65256 0.00000 -0.90296 -1.05960 -0.89413 -1.32883 -1.27830 0.00000 -0.49626 -2.50461 -2.64732 -3.26620 -3.51479 0.00000 octave:19> % A=D*X octave:19> % X = D\A octave:19> % 49*x=49 octave:19> x = (1/49)*49 x = 1.00000 octave:20> format long e octave:21> x x = 1.00000000000000e+00 octave:22> x==1 ans = 0.00000000000000e+00 octave:23> x=(1/3)*3 x = 1.00000000000000e+00 octave:24> x==1 ans = 1.00000000000000e+00 octave:25> 1/49 ans = 2.04081632653061e-02 octave:26> 1/3 ans = 3.33333333333333e-01 octave:27> % 49*x=49 octave:27> x=49/49 x = 1.00000000000000e+00 octave:28> x==1 ans = 1.00000000000000e+00 octave:29> J_J=(speye(5)-D\A) % SI J_J = Columns 1 through 4: 0.00000000000000e+00 -3.77616983459512e+00 -2.11247769243860e+00 -1.26205058942478e-01 -6.68435265434549e-01 0.00000000000000e+00 -5.88409803724241e-01 -9.33408351534935e-01 -8.75023132146121e-01 -6.52556361236922e-01 0.00000000000000e+00 -9.02964625539490e-01 -8.94127506808865e-01 -1.32882596088764e+00 -1.27830486255086e+00 0.00000000000000e+00 -2.50461030165282e+00 -2.64731538052359e+00 -3.26620179860693e+00 -3.51479386829378e+00 Column 5: -3.94262954338153e+00 -1.42085525168871e-01 -1.05959556734569e+00 -4.96261657839181e-01 0.00000000000000e+00 octave:30> format short e octave:31> J_J=(speye(5)-D\A) % SI J_J = 0.0000e+00 -3.7762e+00 -2.1125e+00 -1.2621e-01 -3.9426e+00 -6.6844e-01 0.0000e+00 -5.8841e-01 -9.3341e-01 -1.4209e-01 -8.7502e-01 -6.5256e-01 0.0000e+00 -9.0296e-01 -1.0596e+00 -8.9413e-01 -1.3288e+00 -1.2783e+00 0.0000e+00 -4.9626e-01 -2.5046e+00 -2.6473e+00 -3.2662e+00 -3.5148e+00 0.0000e+00 octave:32> help jacobi 'jacobi' is a function from the file /home/accounts/personale/clrmrc90/aa1819/calcolo_numerico2/jacobi.m [x, its, norm_r] = jacobi (A, b, tol, max_its, x0) Compute the solution of A*x=b through the Jacobi method. tol is a vector of two components with the relative and the absolute tolerances (defaults [1e-6, 1e-6]). max_its is the maximum number of iterations (default 100) and x0 the initial solution (default zero). x is the computed solution, its the number of iterations and norm_r the norm of the final residual. Additional help for built-in functions and operators is available in the online version of the manual. Use the command 'doc ' to search the manual index. Help and information about Octave is also available on the WWW at http://www.octave.org and via the help@octave.org mailing list. octave:33> octave:33> test jacobi warning: Impossibile raggiungere la tolleranza richiesta PASSES 2 out of 2 tests octave:34> A = [2, 1, 0; 1, 2, 1; 0, 1, 2]; octave:35> b = A * ones (3, 1); octave:36> x = jacobi (A, b, 1e-8); octave:37> x x = 1.0000e+00 1.0000e+00 1.0000e+00 octave:38> A = [2, 1, 0; 1, 2, 1; 0, 1, 2]; octave:39> J_J=(speye(3)-diag(diag(A))\A) J_J = 0.0000e+00 -5.0000e-01 0.0000e+00 -5.0000e-01 0.0000e+00 -5.0000e-01 0.0000e+00 -5.0000e-01 0.0000e+00 octave:40> eig(J_J) ans = -7.0711e-01 -4.0332e-17 7.0711e-01 octave:41> A = [1, 1, 1; 1, 2, 1; 0, 1, 1]; octave:42> J_J=(speye(3)-diag(diag(A))\A) J_J = 0.0000e+00 -1.0000e+00 -1.0000e+00 -5.0000e-01 0.0000e+00 -5.0000e-01 0.0000e+00 -1.0000e+00 0.0000e+00 octave:43> eig(J_J) ans = -1.1915e+00 + 0.0000e+00i 5.9574e-01 + 2.5443e-01i 5.9574e-01 - 2.5443e-01i octave:44> ls 2018-10-23.txt 2018-11-06.txt jacobi.m jacobi.m~ octave:45> ls 2018-10-23.txt 2018-11-06.txt gradott.m jacobi.m jacobi.m~ octave:46> A = [1, -1;1, 1]; octave:47> eig(A) ans = 1.0000e+00 + 1.0000e+00i 1.0000e+00 - 1.0000e+00i octave:48> A A = 1.0000e+00 -1.0000e+00 1.0000e+00 1.0000e+00 octave:49> b = A * (1:2)'; octave:50> x0 = [100; 100]; octave:51> [x, its] = gradott (A, b, [], [], x0); warning: Impossibile raggiungere la tolleranza richiesta octave:52> its its = 1.0000e+02 octave:53> x x = 1.0000e+02 1.0000e+02 octave:54> A = [1, 0;1, -1]; octave:55> A A = 1.0000e+00 0.0000e+00 1.0000e+00 -1.0000e+00 octave:56> b = A * (1:2)'; octave:57> x0 = [100; 100]; octave:58> [x, its] = gradott (A, b, [], [], x0); warning: La matrice non sembra definita positiva warning: La matrice non sembra definita positiva warning: La matrice non sembra definita positiva warning: La matrice non sembra definita positiva warning: La matrice non sembra definita positiva octave:59> x x = 1.0000e+00 2.0000e+00 octave:60> A*x ans = 1.0000e+00 -1.0000e+00 octave:61> b b = 1.0000e+00 -1.0000e+00 octave:62> [x, its] = gradott (A, b, [], [], x0); error: La matrice non sembra definita positiva error: called from: error: /home/accounts/personale/clrmrc90/aa1819/calcolo_numerico2/gradott.m at line 33, column 6 octave:62> H=hilb(20); octave:63> cond(H) ans = 2.5703e+18 octave:64> chol(H) error: chol: input matrix must be positive definite octave:64> b=H*(1:20)'; octave:65> H\b warning: matrix singular to machine precision, rcond = 1.99525e-19 ans = 1.0000e+00 2.0000e+00 3.0010e+00 3.9858e+00 5.1039e+00 5.5771e+00 7.9549e+00 6.9938e+00 8.9560e+00 1.0851e+01 1.1093e+01 1.1272e+01 1.2570e+01 1.4475e+01 1.5673e+01 1.6010e+01 1.6265e+01 1.7610e+01 1.9951e+01 1.9660e+01 octave:66> norm(x-(1:20)',inf) error: operator -: nonconformant arguments (op1 is 2x1, op2 is 20x1) error: evaluating argument list element number 1 octave:66> x=H\b x = 1.0000e+00 2.0000e+00 3.0010e+00 3.9858e+00 5.1039e+00 5.5771e+00 7.9549e+00 6.9938e+00 8.9560e+00 1.0851e+01 1.1093e+01 1.1272e+01 1.2570e+01 1.4475e+01 1.5673e+01 1.6010e+01 1.6265e+01 1.7610e+01 1.9951e+01 1.9660e+01 octave:67> norm(x-(1:20)',inf) ans = 1.0062e+00 octave:68> chol(H) error: chol: input matrix must be positive definite octave:68> help gradott 'gradott' is a function from the file /home/accounts/personale/clrmrc90/aa1819/calcolo_numerico2/gradott.m [x, its, norm_r] = gradott (A, b, tol, max_its, x0) Compute the solution of A*x=b through the gradient method. tol is a vector of two components with the relative and the absolute tolerances (defaults [1e-6, 1e-6]). max_its is the maximum number of iterations (default 100) and x0 the initial solution (default zero). x is the computed solution, its the number of iterations and norm_r the norm of the final residual. A has to be symmetric and positive definite. Additional help for built-in functions and operators is available in the online version of the manual. Use the command 'doc ' to search the manual index. Help and information about Octave is also available on the WWW at http://www.octave.org and via the help@octave.org mailing list. octave:69> x=gradott (H, b) warning: Impossibile raggiungere la tolleranza richiesta x = 1.3472e+00 1.5777e+00 1.5408e+00 2.8709e+00 4.6720e+00 6.4992e+00 8.1752e+00 9.6426e+00 1.0896e+01 1.1950e+01 1.2826e+01 1.3549e+01 1.4139e+01 1.4617e+01 1.5000e+01 1.5301e+01 1.5533e+01 1.5708e+01 1.5832e+01 1.5915e+01 octave:70> norm(x-(1:20)',inf) ans = 4.0849e+00 octave:71> x=gradott (H, b, [], 200); warning: Impossibile raggiungere la tolleranza richiesta octave:72> norm(x-(1:20)',inf) ans = 3.0686e+00 octave:73> x=gradott (H, b, [], 1000); warning: Impossibile raggiungere la tolleranza richiesta octave:74> norm(x-(1:20)',inf) ans = 1.7163e+00 octave:75> A = rand (4); octave:76> A=A'*A A = 1.1616e+00 1.0853e+00 7.5562e-01 8.9336e-01 1.0853e+00 1.2950e+00 5.1967e-01 8.5526e-01 7.5562e-01 5.1967e-01 8.2983e-01 8.8349e-01 8.9336e-01 8.5526e-01 8.8349e-01 1.1593e+00 octave:77> b = rand (4, 1); octave:78> x = A \ b; octave:79> E=@(xk) (xk - x)' * A * (xk - x); octave:80> E(rand(4,1)) ans = 3.7221e+01 octave:81> nablaE = @(xk) 2 * xk' * A - 2 * b'; octave:82> nablaE(rand(4,1)) ans = 3.3917e+00 3.9866e+00 1.7462e+00 3.6494e+00 octave:83> v = rand (4, 1); octave:84> xk = rand (4, 1); octave:85> nablaE (xk) * v ans = 3.2984e+00 octave:86> epsilon = 1e-4; octave:87> (E(xk + epsilon * v) - E(xk)) / epsilon ans = 3.2987e+00 octave:88> diary off