hold on %初始参数设置 % i:迭代记数 % e:解精度 % x0:初始解 % x1:迭代一次 ?0:x0点 Jacobi 矩阵 ?1:x1点 Jacobi 矩阵 i=0; e=0; df1=0; x1=0; x0=a; f0=feval(f,x0);d=ones(size(f0'));
df0=Jacbi_div(f,x0); %调用实验<六> Jacbi_div求Jacobi矩阵 while sqrt((f0.^2)*d)>=daltf
x1=x0-f0*inv(df0'); %新迭代得到值,x0为行向量 f1=feval(f,x1); df1=Jacbi_div(f,x1); i=i+1;
if sqrt(x1.^2*d)<1 e=sqrt((x0-x1).^2*d); else
e=sqrt(((x0-x1).^2)*d./(x1.^2))*d; end if e end %满足精度,退出 if i==100 disp('No root find!'); break end %安全设置(最大迭代次数:100) if flag==1 plot([x0 x1],[f0 0],’r’,[x0 x0],[0,f0],’-.r’); %画图 end x0=x1; df0=df1; f0=f1; %准备下次迭代 disp(i) disp([x0' f0']) %输出跌代数据 end 26 if flag==1 X=x1-(a-x1):(a-x1)/50:x1+(a-x1); plot( ,feval(f, ),'k'); %画函数图象 end hold off %输出设置 val=x1; if nargout==2 n=i; end %记时结束 disp('耗时:') disp(toc) 例1: 用 Newton切线法求非线性方程 y?1?e?x?0在0.2 附近的根并画出迭代过程图形. 步骤: 1)新建函数fun1.m function y=fun1(x) %待界解方程 y=1-exp(-x.^2); 2)运算,显示跌代过程如下: >> [val,n]=Newton_root('fun1',0.2 ,1e-6,1e-8) 1 0.0979730645190195 0.00955280068932718 2 0.0487506741811492 0.00237380628825135 3 0.0243463485729836 0.000592569050408831 4 0.012169565781348 0.000148087365290039 5 0.00608433229533475 3.70184142816088e-005 6 0.00304210983785066 9.25438944343604e-006 7 0.00152104788065106 2.31358397884129e-006 8 0.00076052306057252 5.7839515843483e-007 耗时: 0.109999999999999 val = 0.00076052306057252 n = 8 27 2数值分析实验七 第八章 矩阵的特征值计算 1. 实验目的 1) 幂法求矩阵的特征值 2) Jacobi迭代法求矩阵的特征值 2. 实验内容 (1) 幂法求计算矩阵的主特征值及主特征向量 幂法是一种计算矩阵的主特征值及主特征向量的一种跌代法.过程中矩阵不变,适合稀疏矩阵,单有时收敛过慢. function [x0,v,n]=Mifa_eig(A,dalt,x0,p) % P223 幂法求主特征向量及特征值 %[x,v,n]=Mifa_eig(A,dalt,x0,p) % x:输出主特征向量 % v:输出主特征向量 % n:迭代次数 % A :方阵 últ:精度 % x0:初始向量 % p:中心加速 % P223 幂法求主特征向量及特征值 n=size(A); n=n(1); e=1; v=0; m=0; if nargin<4 p=0; end if nargin<3 x0=ones(n,1); end if nargin<2 dalt=1e-6; end for i=1:n 28 A (i,i)=A (i,i)-p; end while e>dalt x=A*x0; t=0; for i=1:n if abs(x(i))>t t=x(i); end end x=x/t; e=max(abs(t-v),max(abs(x-x0))); v=t; x0=x; m=m+1; end v=v+p; if nargout>2 n=m; end 例1 产生一对称矩阵,对不同的原点平移和初值用幂法求计算矩阵的主特征值及主特征向量. >>A=rand(4,4)*10; >>A=A+A'; >>x0=[2 1 3 1]'; >>[x,v,n]=Mifa_eig(A) >>[x,v,n]=Mifa_eig(A,1e-10,x0) >>[x,v,n]=Mifa_eig(A,1e-10,[1 1 1 1]',1.5) >>[x,v,n]=Mifa_eig(A,1e-10,x0,3) (2) Jacobi 方法计算实对称矩阵所有特征值和特征向量 Jacobi 方法通过一组平面旋转 (正交相似变换)将对称矩阵化为对角阵, 得所有特征值和特征向量. 若A?Rn?n为对角阵, 则存在正交阵P,使 PAPT?giag??1,?2,?3, function [v,R]=Jacobi_eigs(A,dalt) %用 Jacobi方法计算实对称矩阵所有特征值 % 和特征向量 % [v,R]=Jacobi_eigs(A,dalt) 29 ,?n??D % v 所有特征值 % R 阵交变换矩阵 % A 对角化矩阵 % dalt 精度 tic if A~=A' error('输入对称方阵!'); end if nargin<2 dalt=1e-6; end n=size(A); n=n(2); e=0; R=eye(n); m=0; for i=1:n for j=i+1:n e=e+A (i,j)^2; end end e=e; while e>=dalt m=m+1; i=1; j=2; for c=1:n for s=c+1:n if abs(A (c,s))>abs(A (i,j)) i=c; j=s; end end end e=e-A (i,j)^2; d=(A (i,i)-A (j,j))/(2*A (i,j)); t=1/(abs(d)+sqrt(d^2+1)); if d>0 c=1/sqrt(1+t^2); s=c*t; 30 else c=-1/sqrt(1+t^2); s=-c*t; end p=A (i,:); q=A (j,:); A (i,i)=p(i)*c^2+q (j)*s^2+2*p(j)*s*c; A (j,j)=p(i)*s^2+q (j)*c^2-2*p(j)*s*c; A (i,j)=(q (j)-p(i))*c*s+p(j)*(c^2-s^2); A (j,i)=A (i,j); for t=1:n if t~=i&t~=j A (i,t)=p(t)*c+q (t)*s; A (t,i)=A (i,t); A (j,t)=q (t)*c-p(t)*s; A (t,j)=A (j,t); end end p=R(:,i); q=R(:,j); for t=1:n R(t,i)=p(t)*c+q (t)*s; R(t,j)=-p(t)*s+q (t)*c; end if m>200 %break; end end v=A (1,1); for i=2:n v=[v A (i,i)]; end toc 31