程序并不显示的计算出H。
Householder变换阵H通用程序: function Householder(A) [n,~]=size(A); e1=zeros(n,1); e1(1)=1; x=A(:,1); w=x-norm(x)*e1; U=w/sqrt((w'*w)); H=eye(n)-2*(U*U'); HA=H*A; disp('HA='); disp(HA); 考虑矩阵
?1??-1A??-2???10?0?24??323?2e??
?2?37?2752??3用你编制的程序计算H使得HA的第一列为?e1的形式,并将HA的结果显示:
HA的结果显示的通用程序:
function Householder(A) [n,~]=size(A); e1=zeros(n,1); e1(1)=1; x=A(:,1);
w=x-norm(x)*e1; U=w/sqrt((w'*w)); H=eye(n)-2*(U*U'); HA=H*A; disp('HA='); disp(HA); 程序运行结果:
>> A=[1 2 3 4;-1 3 sqrt(2) sqrt(3);-2 2 exp(1) pi;-sqrt(10) 2 -3 7;0 2 7 5/2]; >> Householder(A) HA=
4.0000 -2.8311 1.4090 -6.5378 0.0000 1.3896 0.8839 -1.7805 0.0000 -1.2208 1.6576 -3.8836 0.0000 -3.0925 -4.6770 -4.1078 0 2.0000 7.0000 2.5000
7:用jacobi和Gauss-Scidel迭代求解下面的方程组,输出每一步的误差xk?x*;
?5x1?x2?3x3??2???x1?2x2?4x3?1 ??3x?4x?15x?10123?Jacobi迭代通用程序: function jacobi(A,b,x0,ep) if nargin==3 ep=1.0e-6; else if nargin< 3 error return
end end
D=diag(diag(A)); L=-tril(A,-1); U=-triu(A,1); B=D\\(L+U); f=D\\b; y=B*x0+f; n=1;
while norm(y-x0)>=ep fprintf('%d步误差\\n',n) Error=y-x0; disp (Error) x0=y; y=B*x0+f; n=n+1; end
disp('近似解') disp(y)
在命令窗口输入以下命令就可以输出每一步的误差及最优近似解: >> A=[5 -1 -3;-1 2 4;-3 4 15]; >> b=[-2 1 10]'; >> x0=[0 0 0]'; >> ep=10e-8; >> jacobi(A,b,x0,ep)
近似解
-0.081967231207363 -1.803278647380326 1.131147556193974 Gauss-Scidel迭代法通用程序: function Gaussseidel(A,b,x0,ep) if nargin==3 ep=1.0e-6; else if nargin< 3 error return end end
D=diag(diag(A)); L=-tril(A,-1); U=-triu(A,1); B=(D-L)\\U; f=(D-L)\\b; y=B*x0+f; n=1;
while norm(y-x0)>=ep fprintf('%d步误差\\n',n) Error=y-x0; disp (Error) x0=y; y=B*x0+f;
n=n+1; end
disp('x的近似最优解:') disp(y) ;
在命令窗口输入以下命令就可以输出每一步的误差及最优近似解: >> A=[5 -1 -3;-1 2 4;-3 4 15]; >> b=[-2 1 10]'; >> x0=[0 0 0]'; >> ep=10e-8;
>> Gaussseidel(A,b,x0,ep) x的近似最优解: -0.081967202756916 -1.803278585134911 1.131147515484593
8:取不同的初值用Newton迭代法以及弦截法求解x3?2x2?10x?100?0的实根,列表或者画图说明收敛性;
Newton迭代法通用程序:
%使用说明 f为符号表达式,x0为初始值 function Newton(f,x0) syms x; x=x0;
x1=x0-eval(f/diff(f)); while(norm(x1-x0)>0.000001) x=x1;
xk=x1-eval(f/diff(f)); x0=x1;