Notebook source code: https://gist.github.com/MHenderson/7544475
Notebook viewer: http://nbviewer.ipython.org/7544475
%load_ext octavemagic
%octave addpath(genpath("/home/matthew/workspace/resources/N/Numerical Methods/BEARDAHCC/work"))
%%octave
L = [1 0 0 0;2 1 0 0;3 2 1 0;-1 0 2 1];
b = [1 1 1 1]';
L\b
ans = 1 -1 0 2
%%octave
U = [-2 1 1 -1;0 -1 0 1;0 0 1 1;0 0 0 2];
y = [1 -1 0 2]';
x = U\y
x = -0.50000 2.00000 -1.00000 1.00000
%%octave
A = [-2 1 1 -1;-4 1 2 -1;-6 1 4 0;2 -1 1 5];
[L,U] = lufac(A,1)
L = 1 0 0 0 2 1 0 0 3 2 1 0 -1 -0 2 1 U = -2 1 1 -1 0 -1 0 1 0 0 1 1 0 0 0 2
%octave help lufac
`lufac' is a function from the file /home/matthew/workspace/resources/N/Numerical Methods/BEARDAHCC/work/na/lufac.m LUFAC LU factorisations. Example: a=[12 2 3;2 5 1;3 1 1] [l,u]=lufac(a,1); [l,u]=lufac(a,3); [L,U]=LUFAC(A,METHOD,D) returns the LU factors of the square matrix A in matrices L and U. METHOD=1 uses Doolittle's factorisation (diagonals of L=1). METHOD=2 uses Crout's factorisation (diagonals of U=1). METHOD=3 uses Cholesky's factorisation. (Matrix A must be symmetric positive definite.) D (optional) specifies the diagonal elements of the matrix L (METHOD=1) or U (METHOD=2). Example: a=[1 2 4;2 -1 1;3 0 -1] d=[1 2 3] [l,u]=lufac(a,1,d); [l,u]=lufac(a,2,d); Christian C. Beardah 1994 Additional help for built-in functions and operators is available in the on-line version of the manual. Use the command `doc <topic>' 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
A = [-2 1 1 -1;-4 1 2 -1;-6 1 4 0;2 -1 1 5];
b = [1 1 1 1]';
[L,U] = lufac(A, 1)
y = L\b
x = U\y
L = 1 0 0 0 2 1 0 0 3 2 1 0 -1 -0 2 1 U = -2 1 1 -1 0 -1 0 1 0 0 1 1 0 0 0 2 y = 1 -1 0 2 x = -0.50000 2.00000 -1.00000 1.00000
%%octave
newton('newtp5',1,30,1e-3,2,0);
i= 0 p_i= 1.00000 f(p_i)= 144.00000 i= 1 p_i= -0.09091 f(p_i)= 69.48616 i= 2 p_i= -165.99875 f(p_i)= -133690914784.86263 i= 3 p_i= -132.41222 f(p_i)= -43805173332.59692 i= 4 p_i= -105.54626 f(p_i)= -14352720936.39162 i= 5 p_i= -84.05756 f(p_i)= -4702408318.80671 i= 6 p_i= -66.87164 f(p_i)= -1540533854.57490 i= 7 p_i= -53.12917 f(p_i)= -504623928.00724 i= 8 p_i= -42.14293 f(p_i)= -165264972.95593 i= 9 p_i= -33.36351 f(p_i)= -54108499.96578 i= 10 p_i= -26.35174 f(p_i)= -17707369.37320 i= 11 p_i= -20.75676 f(p_i)= -5790872.57505 i= 12 p_i= -16.29839 f(p_i)= -1891825.80747 i= 13 p_i= -12.75308 f(p_i)= -617068.38919 i= 14 p_i= -9.94264 f(p_i)= -200792.65999 i= 15 p_i= -7.72528 f(p_i)= -65098.88986 i= 16 p_i= -5.98851 f(p_i)= -20984.17174 i= 17 p_i= -4.64375 f(p_i)= -6698.62019 i= 18 p_i= -3.62275 f(p_i)= -2099.42838 i= 19 p_i= -2.87611 f(p_i)= -631.86863 i= 20 p_i= -2.37360 f(p_i)= -171.30304 i= 21 p_i= -2.09846 f(p_i)= -34.36219 i= 22 p_i= -2.00905 f(p_i)= -2.87678 i= 23 p_i= -2.00009 f(p_i)= -0.02685 i= 24 p_i= -2.00000 f(p_i)= -0.00000 Converged
%%octave
type newtp5
newtp5 is the user-defined function defined from: /home/matthew/workspace/resources/N/Numerical Methods/BEARDAHCC/work/fnm1/newtp5.m function [y,dydx]=newtp5(x) y=x.^5 - 10*x.^4 + 12*x.^3 + 60*x.^2 + 11*x + 70; dydx=5*x.^4 - 40*x.^3 + 36*x.^2 + 120*x + 11;
%%octave
newton('newtq4',2+2i,30,1e-3,2,0);
i= 1 real(p_i)= 2.00000 imag(p_i)= 2.00000 mod(f(p_i))= 26.83282 i= 2 real(p_i)= 1.94706 imag(p_i)= 1.48824 mod(f(p_i))= 7.67789 i= 3 real(p_i)= 1.95003 imag(p_i)= 1.17388 mod(f(p_i))= 1.90468 i= 4 real(p_i)= 1.97882 imag(p_i)= 1.03021 mod(f(p_i))= 0.31722 i= 5 real(p_i)= 1.99805 imag(p_i)= 1.00049 mod(f(p_i))= 0.01657 i= 6 real(p_i)= 2.00000 imag(p_i)= 0.99999 mod(f(p_i))= 0.00005 Converged
%%octave
A7 = [3 0 1;0 -3 -0;1 0 3]
i = 0
x = [1;1;1;]
i = i + 1, x = A7*x
A7 = 3 0 1 0 -3 -0 1 0 3 i = 0 x = 1 1 1 i = 1 x = 4 -3 4
%%octave
i = i + 1, x = A7*x
i = i + 1, x = A7*x
i = i + 1, x = A7*x
i = i + 1, x = A7*x
i = 2 x = 16 9 16 i = 3 x = 64 -27 64 i = 4 x = 256 81 256 i = 5 x = 1024 -243 1024
%%octave
A7 = [3 0 1;0 -3 -0;1 0 3]
i = 0
x = [1;1;1;]
i = i + 1, y = A7*x, [z,j]=max(abs(y)), delta=y(j), x=y/delta
A7 = 3 0 1 0 -3 -0 1 0 3 i = 0 x = 1 1 1 i = 1 y = 4 -3 4 z = 4 j = 1 delta = 4 x = 1.00000 -0.75000 1.00000
%%octave
i = i + 1, y = A7*x, [z,j]=max(abs(y)), delta=y(j), x=y/delta
i = i + 1, y = A7*x, [z,j]=max(abs(y)), delta=y(j), x=y/delta
i = 2 y = 4.00000 2.25000 4.00000 z = 4 j = 1 delta = 4 x = 1.00000 0.56250 1.00000 i = 3 y = 4.00000 -1.68750 4.00000 z = 4 j = 1 delta = 4 x = 1.00000 -0.42188 1.00000
%%octave
load matrices
x0 = rand(3,1)
[delta,x,n]=scalepm(M3, x0, 200, 0.001);
x0 = 0.12787 0.28723 0.26850
%%octave
eig(M3)
ans = 3.99576 - 0.26780i 4.02022 + 2.48393i -5.01598 - 0.21613i
%%octave
A7 = [3 0 1;0 -3 0;1 0 3]
A7 = 3 0 1 0 -3 0 1 0 3
%%octave
B = inv(A7)
B = 0.37500 0.00000 -0.12500 0.00000 -0.33333 0.00000 -0.12500 0.00000 0.37500
%%octave
i = 0
x = [0;1;2]
i=i+1, y=B*x, [z,j]=max(abs(y)), delta=y(j), x=y/delta
i = 0 x = 0 1 2 i = 1 y = -0.25000 -0.33333 0.75000 z = 0.75000 j = 3 delta = 0.75000 x = -0.33333 -0.44444 1.00000
%%octave
i=i+1, y=B*x, [z,j]=max(abs(y)), delta=y(j), x=y/delta
i = 2 y = -0.25000 0.14815 0.41667 z = 0.41667 j = 3 delta = 0.41667 x = -0.60000 0.35556 1.00000
%%octave
i=i+1, y=B*x, [z,j]=max(abs(y)), delta=y(j), x=y/delta
i = 3 y = -0.35000 -0.11852 0.45000 z = 0.45000 j = 3 delta = 0.45000 x = -0.77778 -0.26337 1.00000
%%octave
i=i+1, y=B*x, [z,j]=max(abs(y)), delta=y(j), x=y/delta
i = 4 y = -0.41667 0.08779 0.47222 z = 0.47222 j = 3 delta = 0.47222 x = -0.88235 0.18591 1.00000
%%octave
i=i+1, y=B*x, [z,j]=max(abs(y)), delta=y(j), x=y/delta
i = 5 y = -0.45588 -0.06197 0.48529 z = 0.48529 j = 3 delta = 0.48529 x = -0.93939 -0.12770 1.00000
%%octave
i=i+1, y=B*x, [z,j]=max(abs(y)), delta=y(j), x=y/delta
i = 6 y = -0.47727 0.04257 0.49242 z = 0.49242 j = 3 delta = 0.49242 x = -0.96923 0.08644 1.00000
%%octave
A7 = [3 0 1;0 -3 0;1 0 3]
i = 0
x = [0;1;2]
i=i+1; y=A7\x, [z,j]=max(abs(y)), delta=y(j), x=y/delta
A7 = 3 0 1 0 -3 0 1 0 3 i = 0 x = 0 1 2 y = -0.25000 -0.33333 0.75000 z = 0.75000 j = 3 delta = 0.75000 x = -0.33333 -0.44444 1.00000
%%octave
i=i+1; y=A7\x, [z,j]=max(abs(y)), delta=y(j), x=y/delta
y = -0.25000 0.14815 0.41667 z = 0.41667 j = 3 delta = 0.41667 x = -0.60000 0.35556 1.00000
%%octave
i=i+1; y=A7\x, [z,j]=max(abs(y)), delta=y(j), x=y/delta
y = -0.35000 -0.11852 0.45000 z = 0.45000 j = 3 delta = 0.45000 x = -0.77778 -0.26337 1.00000
%%octave
i=i+1; y=A7\x, [z,j]=max(abs(y)), delta=y(j), x=y/delta
y = -0.41667 0.08779 0.47222 z = 0.47222 j = 3 delta = 0.47222 x = -0.88235 0.18591 1.00000
%%octave
i=i+1; y=A7\x, [z,j]=max(abs(y)), delta=y(j), x=y/delta
y = -0.45588 -0.06197 0.48529 z = 0.48529 j = 3 delta = 0.48529 x = -0.93939 -0.12770 1.00000
%%octave
i=i+1; y=A7\x, [z,j]=max(abs(y)), delta=y(j), x=y/delta
y = -0.47727 0.04257 0.49242 z = 0.49242 j = 3 delta = 0.49242 x = -0.96923 0.08644 1.00000
%%octave
i=i+1; y=A7\x, [z,j]=max(abs(y)), delta=y(j), x=y/delta
y = -0.48846 -0.02881 0.49615 z = 0.49615 j = 3 delta = 0.49615 x = -0.98450 -0.05807 1.00000
%%octave
i=i+1; y=A7\x, [z,j]=max(abs(y)), delta=y(j), x=y/delta
y = -0.49419 0.01936 0.49806 z = 0.49806 j = 3 delta = 0.49806 x = -0.99222 0.03887 1.00000
%%octave
i=i+1; y=A7\x, [z,j]=max(abs(y)), delta=y(j), x=y/delta
y = -0.49708 -0.01296 0.49903 z = 0.49903 j = 3 delta = 0.49903 x = -0.99610 -0.02596 1.00000
%%octave
A7 =[3 0 1;0 -3 0;1 0 3]
[L,U]=lufac(A7,1)
i=0
x=[0;1;2]
i=i+1; y=U\(L\x), [z,j]=max(abs(y)), delta=y(j), x=y/delta
A7 = 3 0 1 0 -3 0 1 0 3 L = 1.00000 0.00000 0.00000 0.00000 1.00000 0.00000 0.33333 -0.00000 1.00000 U = 3.00000 0.00000 1.00000 0.00000 -3.00000 0.00000 0.00000 0.00000 2.66667 i = 0 x = 0 1 2 y = -0.25000 -0.33333 0.75000 z = 0.75000 j = 3 delta = 0.75000 x = -0.33333 -0.44444 1.00000
%%octave
i=i+1; y=U\(L\x), [z,j]=max(abs(y)), delta=y(j), x=y/delta
y = -0.25000 0.14815 0.41667 z = 0.41667 j = 3 delta = 0.41667 x = -0.60000 0.35556 1.00000
%%octave
i=i+1; y=U\(L\x), [z,j]=max(abs(y)), delta=y(j), x=y/delta
y = -0.35000 -0.11852 0.45000 z = 0.45000 j = 3 delta = 0.45000 x = -0.77778 -0.26337 1.00000
%%octave
i=i+1; y=U\(L\x), [z,j]=max(abs(y)), delta=y(j), x=y/delta
y = -0.41667 0.08779 0.47222 z = 0.47222 j = 3 delta = 0.47222 x = -0.88235 0.18591 1.00000
%%octave
i=i+1; y=U\(L\x), [z,j]=max(abs(y)), delta=y(j), x=y/delta
y = -0.45588 -0.06197 0.48529 z = 0.48529 j = 3 delta = 0.48529 x = -0.93939 -0.12770 1.00000
%%octave
i=i+1; y=U\(L\x), [z,j]=max(abs(y)), delta=y(j), x=y/delta
y = -0.47727 0.04257 0.49242 z = 0.49242 j = 3 delta = 0.49242 x = -0.96923 0.08644 1.00000
%%octave
i=i+1; y=U\(L\x), [z,j]=max(abs(y)), delta=y(j), x=y/delta
y = -0.48846 -0.02881 0.49615 z = 0.49615 j = 3 delta = 0.49615 x = -0.98450 -0.05807 1.00000
%%octave
i=i+1; y=U\(L\x), [z,j]=max(abs(y)), delta=y(j), x=y/delta
y = -0.49419 0.01936 0.49806 z = 0.49806 j = 3 delta = 0.49806 x = -0.99222 0.03887 1.00000
%%octave
A = rand(3)
sigma=1
[V,D]=eig(A)
eig(A)-sigma
[V,D]=eig(A-sigma*eye(3))
A = 0.69417 0.94318 0.01898 0.39730 0.79071 0.68434 0.98099 0.60097 0.95634 sigma = 1 V = -0.41581 + 0.00000i 0.66714 + 0.00000i 0.66714 - 0.00000i -0.55327 + 0.00000i -0.32183 + 0.32725i -0.32183 - 0.32725i -0.72180 + 0.00000i -0.33813 - 0.47950i -0.33813 + 0.47950i D = Diagonal Matrix 1.98210 + 0.00000i 0 0 0 0.22956 + 0.44901i 0 0 0 0.22956 - 0.44901i ans = 0.98210 + 0.00000i -0.77044 + 0.44901i -0.77044 - 0.44901i V = 0.41581 + 0.00000i 0.66714 + 0.00000i 0.66714 - 0.00000i 0.55327 + 0.00000i -0.32183 + 0.32725i -0.32183 - 0.32725i 0.72180 + 0.00000i -0.33813 - 0.47950i -0.33813 + 0.47950i D = Diagonal Matrix 0.98210 + 0.00000i 0 0 0 -0.77044 + 0.44901i 0 0 0 -0.77044 - 0.44901i
%%octave
A7=[3 0 1;0 -3 0;1 0 3]
sigma=-2.5;
M=A7-sigma*eye(3)
i=0
x=rand(3,1)
A7 = 3 0 1 0 -3 0 1 0 3 M = 5.50000 0.00000 1.00000 0.00000 -0.50000 0.00000 1.00000 0.00000 5.50000 i = 0 x = 0.20722 0.87005 0.59929
%%octave
i=i+1, y=M\x, [z,j]=max(abs(y)), delta=y(j), x=y/delta
i = 1 y = 0.01848 -1.74010 0.10560 z = 1.7401 j = 2 delta = -1.7401 x = -0.01062 1.00000 -0.06069
%%octave
i=i+1, y=M\x, [z,j]=max(abs(y)), delta=y(j), x=y/delta
i = 2 y = 0.00008 -2.00000 -0.01105 z = 2 j = 2 delta = -2 x = -0.00004 1.00000 0.00552
%%octave
i=i+1, y=M\x, [z,j]=max(abs(y)), delta=y(j), x=y/delta
i = 3 y = -0.00020 -2.00000 0.00104 z = 2 j = 2 delta = -2 x = 0.00010 1.00000 -0.00052
%%octave
M2=[4 -1 -1 -1;-1 4 -1 -1;-1 -1 4 -1;-1 -1 -1 4]
x0=rand(4,1)
[lambda,x,n]=shiftip(M2,2.9,x0,100,1e-3);
M2 = 4 -1 -1 -1 -1 4 -1 -1 -1 -1 4 -1 -1 -1 -1 4 x0 = 0.82653 0.87984 0.59299 0.29283
%%octave
lambda
lambda = Columns 1 through 8: 0.94009 0.70146 1.15294 0.74514 1.13211 0.78379 1.11327 0.81759 Columns 9 through 16: 1.09647 0.84683 1.08166 0.87192 1.06877 0.89327 1.05766 0.91133 Columns 17 through 24: 1.04815 0.92652 1.04007 0.93923 1.03326 0.94984 1.02754 0.95865 Columns 25 through 32: 1.02276 0.96596 1.01878 0.97200 1.01547 0.97699 1.01273 0.98111 Columns 33 through 40: 1.01047 0.98449 1.00860 0.98728 1.00706 0.98957 1.00579 0.99145 Columns 41 through 48: 1.00475 0.99299 1.00390 0.99426 1.00319 0.99530 1.00262 0.99615 Columns 49 through 56: 1.00214 0.99684 1.00176 0.99742 1.00144 0.99788 1.00118 0.99827 Columns 57 through 64: 1.00097 0.99858 1.00079 0.99884 1.00065 0.99905 1.00053 0.99922 Columns 65 through 70: 1.00043 0.99936 1.00036 0.99948 1.00029 0.99957
%%octave
circle([4 0],2)
%%octave
hold on
plot(3.9958, -0.2678, 'x')
%%octave
M3 = [4 i -1;2 -5 i;1 0 4+2i]
x0=[1;1;1]
[lambda,x,n]=shiftip(M3,-5,x0,10,1e-3)
M3 = 4 + 0i 3 + 0i -1 + 0i 2 + 0i -5 + 0i 3 + 0i 1 + 0i 0 + 0i 4 + 2i x0 = 1 1 1 lambda = -6.14512 + 0.10610i -5.49978 - 0.02798i -5.53065 - 0.01943i -5.52915 - 0.02000i x = Columns 1 through 4: 1.00000 + 0.00000i -0.46829 + 0.01463i -0.30332 - 0.00037i -0.31199 - 0.00004i 1.00000 + 0.00000i 1.00000 - 0.00000i 1.00000 - 0.00000i 1.00000 - 0.00000i 1.00000 + 0.00000i -0.06951 + 0.02561i 0.03562 - 0.00908i 0.03111 - 0.00645i Column 5: -0.31151 - 0.00004i 1.00000 - 0.00000i 0.03129 - 0.00664i n = 4