<P>function I=R_quadx(s,t)</P>
<P>%输入a,e,mmax<BR>a=1</P>
<P>%离散化<BR>t=0:0.1:1;<BR>n=length(t);<BR>d=(t(n)-t(1))/(n-1)<BR>s=d/2:d:(d/2)+((n-2)*d);<BR>t=s;<BR>n=length(t);</P>
<P>for i=1:n<BR> for j=1:n<BR> A(i,j)=d*exp(-a*abs(s(i)-t(j)));<BR> end<BR> b(i)=(2*s(i)/a)-(exp((-a)*(1-s(i)))/a)+((1/(a*a))*(exp(-a*s(i))-exp(-a*(1-s(i)))));<BR>end</P>
<P>for i=1:n<BR> G(i)=A(i,1)*t(1);<BR> for j=1:(n-1)<BR> G(i)=G(i)+A(i,j+1)*t(j+1);<BR> end<BR>end</P>
<P>[m,n]=size(A)<BR>nb=length(b)<BR>ng=length(G)</P>
<P>if n~=m<BR> error('The rows and columns of matrix A must be equal!');<BR> return;<BR>end<BR>if ng~=n<BR> error('The columns of A must be equal the length of G!');<BR> return;<BR>end<BR>if nb~=n<BR> error('The columns of A must be equal the length of b!');<BR> return;<BR>end</P>
<P>A=A<BR>t=t'<BR>b=b'<BR>G=G'<BR>s=s'</P>
<P>%输入x0,m0<BR>x0=t<BR>m0=20</P>
<P>%Arnoldi过程,并输出Hm<BR>r0=b-A*x0<BR>v=r0/norm(r0)<BR>p(:,1)=v;</P>
<P>for j=1:m<BR> f=A*v;<BR> for i=1:j<BR> H(i,j)=dot(f,v);<BR> f=f-H(i,j)*v;<BR> end<BR> H(j+1,j)=norm(f);<BR> v=f/H(j+1,j);<BR> p(:,j+1)=v;<BR>end<BR>H=H<BR>p=p</P>
<P>%对Hm作SVD分解,并且输出{c,v,u}<BR>[U,C,V]=svd(H,0)<BR> </P>
<P>%求解xk及rk=A*xk-b,再进一步求出-log||xk||和-log||rk||,并作图<BR>U(n+1,:)=[];<BR>be=G-b<BR>for k=1:m-1<BR> <BR> for j=1:n<BR> xe(:,j)=dot(be,U(:,j))*(1/C(j,j))*V(:,j);<BR> end<BR> for j=1:n<BR> re(:,j)=dot(be,U(:,j))*V(:,j);<BR> end<BR> xk=xe(:,1);<BR> for j=2:k<BR> xk=xk+xe(:,j);<BR> end<BR> rk=re(:,k+1);<BR> if k+2<=n<BR> for j=k+2:n<BR> rk=rk+re(:,j);<BR> end<BR> e1(k)=(norm(xk))^(-1);<BR> e2(k)=(norm(rk))^(-1);<BR> end<BR>end<BR>loglog(e1,e2)<BR>xe=xe<BR>re=re<BR>xk=xk<BR>rk=rk</P> |