matlab 矩阵和数值分析,浅谈对矩阵理论数值分析及matlab应用的理解 7

浅谈对矩阵理论数值分析及matlab应用的理解 7

附录1 Ax=b 1、 利用高斯消元法 clear at rat A=randn(5) b=randn(5,1) [m,n]=size(A); for i=1:(m-1) numb=int2str(i); for j=(i+1):m A(j,:)=A(j,:)-A(i,:)*A(j,i)/A(i,i); end A end x(m)=A(m,n)/A(m,m); for i=(m-1):-1:1 x(i)=(A(i,n)-A(i,i+1:m)*x(i+1:m) )/A(i,i); end x 2、利用高斯列主元消元法 clear at rat A=randn(5) b=randn(5,1) [m,n]=size(A); for i=1:(m-1) numb=int2str(i); disp([ 第 ,numb, 次选列主元后的增广矩阵 ]) temp=max(abs(A(i:m,i))); [a,b]=find(abs(A(i:m,i))==temp); tempo=A(a(1)+i-1,:); A(a(1)+i-1,:)=A(i,:); A(i,:)=tempo disp([ 第 ,numb, 次消元后的增广矩阵 ]) for j=(i+1):m A(j,:)=A(j,:)-A(i,:)*A(j,i)/A(i,i); end A end %回代过程 disp( 回代求解 ) x(m)=A(m,n)/A(m,m); for i=(m-1):-1:1 x(i)=(A(i,n)-A(i,i+1:m)*x(i+1:m) )/A(i,i); end x 3、LU分解 A = randn(5); [L,U,P] = lu(A) L = P\L 附录2 n =5; A = randn(n); for i = 1:n A(i,i) = sum(abs(A(i,:))); end b = randn(n, 1); D = diag(diag(A)); L=tril(A,-1) U=triu(A,1) M = D\(L+U); f = D\b; Rho = max(abs(eig(M))) R = 1e-07; switch sign(1 - Rho) case -1 disp( The Jacobian is not applicable ) otherwise x(:,1) = normrnd(0, 9, n, 1); k = 1; while k =R k = k+1; else X = x(:, k+1) IterN = k break end end end Y = A\b ER_1 = norm(Y - X, 1) ER_2 = norm(Y - X, 2) ER_inf = norm(Y - X, inf) 改变谱半径的大小 看收敛速度的变化 n =5; A = randn(n) for i = 1:n A(i,i) = sum(abs(A(i,:))); end b = randn(n, 1) D = diag(diag(A)); L=tril(A,-1); U=triu(A,1); M = D\(L+U); f = D\b; Rho = max(abs(eig(M)))/2 R = 1e-07; switch sign(1 - Rho) case -1 disp( The Jacobian is not applicable ) otherwise x(:,1) = normrnd(0, 9, n, 1); k = 1; while k =R k = k+1; else X = x(:, k+1) IterN = k break end end end Y = A\b; ER_1 = norm(Y - X, 1); ER_2 = norm(Y - X, 2); ER_inf = norm(Y - X, inf) ; n =5; A = randn(n) for i = 1:n A(i,i) = sum(abs(A(i,:))); end b = randn(n, 1) D = diag(diag(A)); L=tril(A,-1); U=triu(A,1); M = D\(L+U); f = D\b; Rho = max(abs(eig(M)))/3 R = 1e-07; switch sign(1 - Rho) case -1 disp( The Jacobian is not applicable ) otherwise x(:,1) = normrnd(0, 9, n, 1); k = 1; while k =R k = k+1; else X = x(:, k+1) IterN = k break end end end Y = A\b; ER_1 = norm(Y - X, 1); ER_2 = norm(Y - X, 2); ER_inf = norm(Y - X, inf) ; n = 5; C = normrnd(1, 4, n); [A R] = qr(C); A = A + 3*eye(n); CondA = cond(A, 2) X = ones(n,1); b = A*X Y = A\b D=diag(diag(A)); L=tril(A,-1); U=triu(A,1); for k = 1:100 W(k) = 0.2*k; Gs(:,:, k) = (D - W(k)*L)\((1 - W(k))*D + W(k)*U); Rho(k) = max(abs(eig(Gs(:,:,k)))); end [MinRho,k] = min(Rho) IterN = k w = W(K) GS = (D - w*L)\((1 - w)*D + w*U); FS = (D - w*L)\(w*b); r = 1e-06; 附录3 clear all x = -3:0.1:3; Fun = @(x)- 2*x.^6+10*x.^4 - 2*x.^3 +5; fun = Fun(x); plot(x, fun, x, 0*x, LineWidth , 3) hold on title( Solve to algebraic equation , fontsize , 24) %Real Roots are in [-1, -0.5], [0.5, 1] and [2, 2.5], Respectively a = 1; b = 3; X = (a + b)/2; Beta = (sqrt(5) - 1)/2; x1 = a + (1 - Beta)*(b - a); x2 = a + Beta*(b - a); fa = Fun(a); fb = Fun(b); fx = Fun(X); f1 = Fun(x1); f2 = Fun(x2); xr = 1e-05; fr = 1e-05; N = 0; while abs(fx) >= fr N = N +1 if abs(b-a) < xr X = (a+b)/2 break elseif f1 == 0 X = x1 break elseif f2 == 0 X = x2 break elseif fa*f1 < 0 b = x1; X = (a + b)/2; x1 = a + (1 - Beta)*(b - a); x2 = a + Beta*(b - a); f1 = Fun(x1); f2 = Fun(x2); fx = Fun(X); elseif fb*f2 < 0 a = x2; X = (a + b)/2; x1 = a + (1 - Beta)*(b - a); x2 = a + Beta*(b - a); f1 = Fun(x1); f2 = Fun(x2); fx = Fun(X); else a = x1; b = x2; X = (a + b)/2; x1 = a + (1 - Beta)*(b - a); x2 = a + Beta*(b - a); f1 = Fun(x1); f2 = Fun(x2); fx = Fun(X); end end plot(X, Fun(X), r. , markersize , 30) at long x0 = -3 Fsolve_x = fsolve( @(x)- 2*x.^6+10*x.^4 - 2*x.^3 +5, x0) Fzero_x = fzero( @(x)- 2*x.^6+10*x.^4 - 2*x.^3 +5, x0) 附录4 function F = Experiment_03_Fun(x) F = [x(1).^2+2*x(1)*x(2)-x(2).^2-1; x(1).^2-0.94*x(1)*x(2)-x(2).^2-9]; % Solve the Equations by Newton Iteration F = @(x, y) [x.^2+b*x*y-y^2-1; x.^2+c*x*y+y^2-9]; Fun = [F1; F2]; Fx = diff(Fun, s); %求偏导 Fy = diff(Fun, t); %求偏导 DF = [Fx Fy]; %生成雅可比矩阵 x(:,1) = [2; -3]; %定义初始值 n = 50; %最大迭代数 Rx = 1e-8; %偏差 k = 1; while k