octave:2> A=rand(4) A = 0.651901 0.792023 0.872030 0.384807 0.709346 0.819125 0.673448 0.241153 0.535944 0.150197 0.617478 0.866001 0.439766 0.292357 0.982676 0.021815 octave:3> b=rand(4,1) b = 0.949193 0.437645 0.475117 0.082355 octave:4> x=A\b x = -3.17165 2.14165 0.83167 1.54704 octave:5> A*x ans = 0.949193 0.437645 0.475117 0.082355 octave:6> b b = 0.949193 0.437645 0.475117 0.082355 octave:7> b=A*ones(4,1) b = 2.7008 2.4431 2.1696 1.7366 octave:8> A\b ans = 1.00000 1.00000 1.00000 1.00000 octave:9> [L,U]=lu(A) L = 0.91902 -0.08371 0.50895 1.00000 1.00000 0.00000 0.00000 0.00000 0.75555 1.00000 0.00000 0.00000 0.61996 0.45972 1.00000 0.00000 U = 0.70935 0.81912 0.67345 0.24115 0.00000 -0.46869 0.10866 0.68380 0.00000 0.00000 0.51521 -0.44205 0.00000 0.00000 0.00000 0.44540 octave:10> U\(L\b)) parse error: syntax error >>> U\(L\b)) ^ octave:10> U\(L\b) ans = 1.00000 1.00000 1.00000 1.00000 octave:11> [L,U,P]=lu(A) L = 1.00000 0.00000 0.00000 0.00000 0.75555 1.00000 0.00000 0.00000 0.61996 0.45972 1.00000 0.00000 0.91902 -0.08371 0.50895 1.00000 U = 0.70935 0.81912 0.67345 0.24115 0.00000 -0.46869 0.10866 0.68380 0.00000 0.00000 0.51521 -0.44205 0.00000 0.00000 0.00000 0.44540 P = Permutation Matrix 0 1 0 0 0 0 1 0 0 0 0 1 1 0 0 0 octave:12> octave:12> U\(L\(P*b)) ans = 1.00000 1.00000 1.00000 1.00000 octave:13> help lu 'lu' is a built-in function from the file libinterp/corefcn/lu.cc -- Built-in Function: [L, U] = lu (A) -- Built-in Function: [L, U, P] = lu (A) -- Built-in Function: [L, U, P, Q] = lu (S) -- Built-in Function: [L, U, P, Q, R] = lu (S) -- Built-in Function: [...] = lu (S, THRES) -- Built-in Function: Y = lu (...) -- Built-in Function: [...] = lu (..., "vector") Compute the LU decomposition of A. If A is full subroutines from LAPACK are used and if A is sparse then UMFPACK is used. The result is returned in a permuted form, according to the optional return value P. For example, given the matrix 'a = [1, 2; 3, 4]', [l, u, p] = lu (A) returns l = 1.00000 0.00000 0.33333 1.00000 u = 3.00000 4.00000 0.00000 0.66667 p = 0 1 1 0 The matrix is not required to be square. When called with two or three output arguments and a spare input matrix, 'lu' does not attempt to perform sparsity preserving column permutations. Called with a fourth output argument, the sparsity preserving column transformation Q is returned, such that 'P * A * Q = L * U'. Called with a fifth output argument and a sparse input matrix, 'lu' attempts to use a scaling factor R on the input matrix such that 'P * (R \ A) * Q = L * U'. This typically leads to a sparser and more stable factorization. An additional input argument THRES, that defines the pivoting threshold can be given. THRES can be a scalar, in which case it defines the UMFPACK pivoting tolerance for both symmetric and unsymmetric cases. If THRES is a 2-element vector, then the first element defines the pivoting tolerance for the unsymmetric UMFPACK pivoting strategy and the second for the symmetric strategy. By default, the values defined by 'spparms' are used ([0.1, 0.001]). Given the string argument "vector", 'lu' returns the values of P and Q as vector values, such that for full matrix, 'A (P,:) = L * U', and 'R(P,:) * A (:, Q) = L * U'. With two output arguments, returns the permuted forms of the upper and lower triangular matrices, such that 'A = L * U'. With one output argument Y, then the matrix returned by the LAPACK routines is returned. If the input matrix is sparse then the matrix L is embedded into U to give a return value similar to the full case. For both full and sparse matrices, 'lu' loses the permutation information. See also: luupdate, chol, hess, qr, qz, schur, svd. 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:14> [L,U,P]=lu(A,'vector') L = 1.00000 0.00000 0.00000 0.00000 0.75555 1.00000 0.00000 0.00000 0.61996 0.45972 1.00000 0.00000 0.91902 -0.08371 0.50895 1.00000 U = 0.70935 0.81912 0.67345 0.24115 0.00000 -0.46869 0.10866 0.68380 0.00000 0.00000 0.51521 -0.44205 0.00000 0.00000 0.00000 0.44540 P = 2 3 4 1 octave:15> [L,U,P]=lu(A) L = 1.00000 0.00000 0.00000 0.00000 0.75555 1.00000 0.00000 0.00000 0.61996 0.45972 1.00000 0.00000 0.91902 -0.08371 0.50895 1.00000 U = 0.70935 0.81912 0.67345 0.24115 0.00000 -0.46869 0.10866 0.68380 0.00000 0.00000 0.51521 -0.44205 0.00000 0.00000 0.00000 0.44540 P = Permutation Matrix 0 1 0 0 0 0 1 0 0 0 0 1 1 0 0 0 octave:16> [L,U,p]=lu(A,'vector') L = 1.00000 0.00000 0.00000 0.00000 0.75555 1.00000 0.00000 0.00000 0.61996 0.45972 1.00000 0.00000 0.91902 -0.08371 0.50895 1.00000 U = 0.70935 0.81912 0.67345 0.24115 0.00000 -0.46869 0.10866 0.68380 0.00000 0.00000 0.51521 -0.44205 0.00000 0.00000 0.00000 0.44540 p = 2 3 4 1 octave:17> b b = 2.7008 2.4431 2.1696 1.7366 octave:18> P*b ans = 2.4431 2.1696 1.7366 2.7008 octave:19> b(p) ans = 2.4431 2.1696 1.7366 2.7008 octave:20> b b = 2.7008 2.4431 2.1696 1.7366 octave:21> b(3) ans = 2.1696 octave:22> b([1,3]) ans = 2.7008 2.1696 octave:23> b([3,1]) ans = 2.1696 2.7008 octave:24> p p = 2 3 4 1 octave:25> b(p) ans = 2.4431 2.1696 1.7366 2.7008 octave:26> U\(L\b(p)) ans = 1.00000 1.00000 1.00000 1.00000 octave:27> whos Variables in the current scope: Attr Name Size Bytes Class ==== ==== ==== ===== ===== A 4x4 128 double L 4x4 128 double P 4x4 16 double U 4x4 128 double ans 4x1 32 double b 4x1 32 double p 4x1 32 double x 4x1 32 double Total is 80 elements using 528 bytes octave:28> L*U ans = 0.709346 0.819125 0.673448 0.241153 0.535944 0.150197 0.617478 0.866001 0.439766 0.292357 0.982676 0.021815 0.651901 0.792023 0.872030 0.384807 octave:29> P*A ans = 0.709346 0.819125 0.673448 0.241153 0.535944 0.150197 0.617478 0.866001 0.439766 0.292357 0.982676 0.021815 0.651901 0.792023 0.872030 0.384807 octave:30> A A = 0.651901 0.792023 0.872030 0.384807 0.709346 0.819125 0.673448 0.241153 0.535944 0.150197 0.617478 0.866001 0.439766 0.292357 0.982676 0.021815 octave:31> A=A*A' A = 1.96079 1.79125 1.34004 1.38356 1.79125 1.68582 1.12788 1.21847 1.34004 1.12788 1.44103 0.90527 1.38356 1.21847 0.90527 1.24500 octave:32> R=chol(A) R = 1.40028 1.27921 0.95698 0.98806 0.00000 0.22237 -0.43308 -0.20445 0.00000 0.00000 0.58108 -0.22169 0.00000 0.00000 0.00000 0.42166 octave:33> R'*R ans = 1.96079 1.79125 1.34004 1.38356 1.79125 1.68582 1.12788 1.21847 1.34004 1.12788 1.44103 0.90527 1.38356 1.21847 0.90527 1.24500 octave:34> A A = 1.96079 1.79125 1.34004 1.38356 1.79125 1.68582 1.12788 1.21847 1.34004 1.12788 1.44103 0.90527 1.38356 1.21847 0.90527 1.24500 octave:35> b=A*ones(4,1) b = 6.4756 5.8234 4.8142 4.7523 octave:36> R\(R'\b) ans = 1.00000 1.00000 1.00000 1.00000 octave:37> ls 2017-10-10.txt newton.m newton.m~ octave:38> f=@(x) [x(1)^2+x(2)^2-1;sin(pi*x(1)/2)+x(2)^3]; octave:39> Jf=@(x) [2*x(1),2*x(2);pi/2*cos(pi*x(1)/2),3*x(2)^2]; octave:40> x0=[1;1]; octave:41> tol=1e-10; octave:42> x=newton(f,Jf,x0,tol) x = 0.47610 -0.87939 octave:43> f(x) ans = 8.7594e-12 -2.0567e-11 octave:44> x=newton(f,Jf,x0,tol) ans = 2.4631 ans = 1.2792 ans = 0.51071 ans = 0.26900 ans = 0.042789 ans = 0.0019066 ans = 2.9596e-06 ans = 7.9440e-12 x = 0.47610 -0.87939 octave:45> Jf=@(x) [2*x(1),2*x(2);pi/2*cos(pi*x(1)/2),3*x(2)]; octave:46> x=newton(f,Jf,x0,tol) ans = 0.97522 ans = 0.43877 ans = 1.3484 ans = 0.70823 ans = 0.84673 ans = 0.65438 ans = 1.5234 ans = 1.1133 ans = 4.7996 ans = 2.3620 ans = 1.0883 ans = 0.90994 ans = 0.68490 ans = 2.1280 ans = 1.4846 ans = 0.82483 ans = 2.0240 ans = 1.3343 ans = 0.76439 ans = 2.4731 ans = 2.1579 ans = 1.0517 ans = 3.1975 ans = 4.5502 ans = 7.2681 ans = 33.444 ans = 27.650 ans = 551.16 ans = 2.2296e+04 ans = 3.1242e+08 ans = 6.7198e+15 ans = 2.8375e+31 ans = 5.5432e+61 ans = 1.9308e+123 ans = Inf warning: matrix singular to machine precision, rcond = nan error: newton: exception encountered in Fortran subroutine dgelsd_ error: called from: error: /home/accounts/personale/clrmrc90/aa1718/equazioni_differenziali/newton.m at line 6, column 11 octave:46> Jf=@(x) [2*x(1),2*x(2);pi/2*cos(pi*x(1)/2),x(2)^2]; octave:47> x=newton(f,Jf,x0,tol) ans = 1.1707 ans = 9.2452 ans = 12.426 ans = 185.66 ans = 2038.3 ans = 2.9238e+05 ans = 4.2354e+07 ans = 3.4037e+14 warning: matrix singular to machine precision, rcond = 3.65633e-22 ans = 3.4037e+14 ans = 3.6738e+09 ans = 6.3726e+11 warning: matrix singular to machine precision, rcond = 0 ans = 3.1863e+11 warning: matrix singular to machine precision, rcond = 0 ans = 1.5932e+11 warning: matrix singular to machine precision, rcond = 0 ans = 7.9658e+10 warning: matrix singular to machine precision, rcond = 0 ans = 3.9829e+10 warning: matrix singular to machine precision, rcond = 0 ans = 1.9914e+10 warning: matrix singular to machine precision, rcond = 0 ans = 9.9572e+09 warning: matrix singular to machine precision, rcond = 0 ans = 4.9786e+09 warning: matrix singular to machine precision, rcond = 0 ans = 2.4893e+09 warning: matrix singular to machine precision, rcond = 0 ans = 1.2446e+09 warning: matrix singular to machine precision, rcond = 0 ans = 6.2232e+08 warning: matrix singular to machine precision, rcond = 0 ans = 3.1116e+08 warning: matrix singular to machine precision, rcond = 0 ans = 1.5558e+08 warning: matrix singular to machine precision, rcond = 0 ans = 7.7791e+07 warning: matrix singular to machine precision, rcond = 0 ans = 3.8895e+07 warning: matrix singular to machine precision, rcond = 0 ans = 1.9448e+07 warning: matrix singular to machine precision, rcond = 0 ans = 9.7238e+06 warning: matrix singular to machine precision, rcond = 0 ans = 4.8619e+06 warning: matrix singular to machine precision, rcond = 0 ans = 2.4310e+06 warning: matrix singular to machine precision, rcond = 0 ans = 1.2155e+06 warning: matrix singular to machine precision, rcond = 0 ans = 6.0774e+05 warning: matrix singular to machine precision, rcond = 0 ans = 3.0387e+05 warning: matrix singular to machine precision, rcond = 0 ans = 1.5193e+05 warning: matrix singular to machine precision, rcond = 0 ans = 7.5967e+04 warning: matrix singular to machine precision, rcond = 0 ans = 3.7984e+04 warning: matrix singular to machine precision, rcond = 0 ans = 1.8992e+04 warning: matrix singular to machine precision, rcond = 0 ans = 9495.9 warning: matrix singular to machine precision, rcond = 0 ans = 4748.0 warning: matrix singular to machine precision, rcond = 0 ans = 2374.0 warning: matrix singular to machine precision, rcond = 0 ans = 1187.0 warning: matrix singular to machine precision, rcond = 0 ans = 593.49 warning: matrix singular to machine precision, rcond = 0 ans = 296.75 warning: matrix singular to machine precision, rcond = 0 ans = 148.37 warning: matrix singular to machine precision, rcond = 0 ans = 74.183 warning: matrix singular to machine precision, rcond = 0 ans = 37.086 warning: matrix singular to machine precision, rcond = 0 ans = 18.540 warning: matrix singular to machine precision, rcond = 0 ans = 9.2507 warning: matrix singular to machine precision, rcond = 0 ans = 4.5953 warning: matrix singular to machine precision, rcond = 0 ans = 2.2503 warning: matrix singular to machine precision, rcond = 0 ans = 1.0095 warning: matrix singular to machine precision, rcond = 0 ans = 0.26384 warning: matrix singular to machine precision, rcond = 0 ans = 0.097749 warning: matrix singular to machine precision, rcond = 0 ans = 0.044679 warning: matrix singular to machine precision, rcond = 0 ans = 0.023222 warning: matrix singular to machine precision, rcond = 0 ans = 0.012963 warning: matrix singular to machine precision, rcond = 0 ans = 0.0075354 warning: matrix singular to machine precision, rcond = 0 ans = 0.0044848 warning: matrix singular to machine precision, rcond = 0 ans = 0.0027069 warning: matrix singular to machine precision, rcond = 0 ans = 0.0016476 warning: matrix singular to machine precision, rcond = 0 ans = 0.0010081 warning: matrix singular to machine precision, rcond = 0 ans = 6.1872e-04 warning: matrix singular to machine precision, rcond = 0 ans = 3.8048e-04 warning: matrix singular to machine precision, rcond = 0 ans = 2.3425e-04 warning: matrix singular to machine precision, rcond = 0 ans = 1.4433e-04 warning: matrix singular to machine precision, rcond = 0 ans = 8.8966e-05 warning: matrix singular to machine precision, rcond = 0 ans = 5.4854e-05 warning: matrix singular to machine precision, rcond = 0 ans = 3.3827e-05 warning: matrix singular to machine precision, rcond = 0 ans = 2.0863e-05 warning: matrix singular to machine precision, rcond = 0 ans = 1.2868e-05 warning: matrix singular to machine precision, rcond = 0 ans = 7.9371e-06 warning: matrix singular to machine precision, rcond = 0 ans = 4.8958e-06 warning: matrix singular to machine precision, rcond = 0 ans = 3.0199e-06 warning: matrix singular to machine precision, rcond = 0 ans = 1.8628e-06 warning: matrix singular to machine precision, rcond = 0 ans = 1.1491e-06 warning: matrix singular to machine precision, rcond = 0 ans = 7.0879e-07 warning: matrix singular to machine precision, rcond = 0 ans = 4.3722e-07 warning: matrix singular to machine precision, rcond = 0 ans = 2.6970e-07 warning: matrix singular to machine precision, rcond = 0 ans = 1.6636e-07 warning: matrix singular to machine precision, rcond = 0 ans = 1.0262e-07 warning: matrix singular to machine precision, rcond = 0 ans = 6.3301e-08 warning: matrix singular to machine precision, rcond = 0 ans = 3.9048e-08 warning: matrix singular to machine precision, rcond = 0 ans = 2.4086e-08 warning: matrix singular to machine precision, rcond = 0 ans = 1.4858e-08 warning: matrix singular to machine precision, rcond = 0 ans = 9.1650e-09 warning: matrix singular to machine precision, rcond = 0 ans = 5.6534e-09 warning: matrix singular to machine precision, rcond = 0 ans = 3.4873e-09 warning: matrix singular to machine precision, rcond = 0 ans = 2.1512e-09 warning: matrix singular to machine precision, rcond = 0 ans = 1.3269e-09 warning: matrix singular to machine precision, rcond = 0 ans = 8.1852e-10 warning: matrix singular to machine precision, rcond = 0 ans = 5.0491e-10 warning: matrix singular to machine precision, rcond = 0 ans = 3.1145e-10 warning: matrix singular to machine precision, rcond = 0 ans = 1.9212e-10 warning: matrix singular to machine precision, rcond = 0 ans = 1.1851e-10 warning: matrix singular to machine precision, rcond = 0 ans = 7.3102e-11 x = 1.00000 0.00000 octave:48> Jf=@(x) [x(1),2*x(2);pi/2*cos(pi*x(1)/2),3*x(2)^2]; octave:49> x=newton(f,Jf,x0,tol) ans = 1.9902 ans = 2.5596 ans = 0.71634 ans = 26.202 ans = 25.480 ans = 11.976 ans = 11.347 ans = 13.098 ans = 12.653 ans = 0.49225 ans = 0.14962 ans = 0.046678 ans = 0.018030 ans = 0.0063414 ans = 0.0022530 ans = 7.9340e-04 ans = 2.8020e-04 ans = 9.8853e-05 ans = 3.4887e-05 ans = 1.2311e-05 ans = 4.3444e-06 ans = 1.5331e-06 ans = 5.4100e-07 ans = 1.9091e-07 ans = 6.7370e-08 ans = 2.3774e-08 ans = 8.3895e-09 ans = 2.9605e-09 ans = 1.0447e-09 ans = 3.6867e-10 ans = 1.3010e-10 ans = 4.5910e-11 x = 0.47610 -0.87939 octave:50> Jf=@(x) [x(1),2*x(2);pi/2*cos(pi*x(1)/2),3*x(2)^2]; octave:51> x=newton(f,Jf,x0,tol) ans = 1.9902 ans = 2.5596 ans = 0.71634 ans = 26.202 ans = 25.480 ans = 11.976 ans = 11.347 ans = 13.098 ans = 12.653 ans = 0.49225 ans = 0.14962 ans = 0.046678 ans = 0.018030 ans = 0.0063414 ans = 0.0022530 ans = 7.9340e-04 ans = 2.8020e-04 ans = 9.8853e-05 ans = 3.4887e-05 ans = 1.2311e-05 ans = 4.3444e-06 ans = 1.5331e-06 ans = 5.4100e-07 ans = 1.9091e-07 ans = 6.7370e-08 ans = 2.3774e-08 ans = 8.3895e-09 ans = 2.9605e-09 ans = 1.0447e-09 ans = 3.6867e-10 ans = 1.3010e-10 ans = 4.5910e-11 x = 0.47610 -0.87939 octave:52> x=newton(f,Jf,x0,tol) x = 0.47610 -0.87939 octave:53> help newton error: help: 'newton' is not documented octave:53> help newton error: help: 'newton' is not documented octave:53> help newton 'newton' is a function from the file /home/accounts/personale/clrmrc90/aa1718/equazioni_differenziali/newton.m x = newton(f,Jf,x0,tol) tol = [relative_tol,absolute_tol] 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:54> x=newton(f,Jf,x0,tol) error: newton: A(I): index out of bounds; value 2 out of bound 1 error: called from: error: /home/accounts/personale/clrmrc90/aa1718/equazioni_differenziali/newton.m at line 8, column 1 octave:54> tol=[1e-10,1e-10] tol = 1.0000e-10 1.0000e-10 octave:55> x=newton(f,Jf,x0,tol) x = 0.47610 -0.87939 octave:56> diary off