《数值分析与算法》例题matlab源程序.doc

数值分析与算法喻文健，清华大学出版社，例题matlab源程序%P33 example 2-2a=1.0;b=1.5; while b-aeps x=a+(b-a)/2 if f_ex22(a)*f_ex22(x)0 a=x; else b=x; endend %P36 example 2-4x=1.5for i=1:10 x=(x+2)(1/4)end %P42 example 2-7x=1.5for i=1:10 x=(3*x4+2)/(4*x3-1)end%P54 example2-9x01=1;x02=2;x0=x01;x02; for i=1:8x01=x0(1,1);x02=x0(2,1); f01=x01+2*x02-2;f02=x012+4*x022-4; J01=1,2;J02=2*x01,8*x02;s0=inv(J01;J02)*f01;f02 x1=x0-s0 x0=x1;endP94 example 3-14A=5 -1 -1;-1 3 -1;-1 -1 5;A(:,1)=A(:,1)./sqrt(A(1,1);A(2,2)=sqrt(A(2,2)-A(2,1)2);A(3,2)=(A(3,2)-A(3,1)2)/A(2,2);A(3,3)=sqrt(A(3,3)-A(3,1)2-A(3,2)2);% P111, 6n=12;a=hilb(n);x=ones(n,1);b=a*x;b=b+10(-7)*ones(n,1);L=chol(a,lower); %Cholesky y=ones(n,1);y(1)=b(1)/L(1,1);for i=2:n y(i)=(b(i)-sum(L(i,1:i-1)*y(1:i-1,1)/L(i,i);% L lowerend M=L;x(n)=y(n)/M(n,n);for i=n-1:-1:1 x(i)=(y(i)-sum(M(i,i+1:n)*y(i+1:n,1)/M(i,i);% L upperend r=b-a*x;delt_x=x-ones(n,1);norm_r=norm(r,inf) % norm of rnorm_delt_x=norm(delt_x,inf) % norm of delt_x% P121 example 4-3A=10 3 1;2 -10 3;1 3 10;b=14 -5 14;x0=0 0 0; % JacobiAA=0 -1/10*3 -1/10;1/5 0 3/10;-1/10 -3/10 0;bb=14/10 5/10 14/10;for i=1:9 x(:,i)=AA*x0+bb; x0=x(:,i);endx %G-Sb=14 -5 14;x0=0 0 0;D=10 0 0;0 -10 0;0 0 10;L=0 0 0;-2 0 0;-1 -3 0;U=0 -3 -1;0 0 -3; 0 0 0;for i=1:9 xx(:,i)=inv(D-L)*U*x0+inv(D-L)*b; x0=xx(:,i);endxx %SORb=14 -5 14;x0=0 0 0;D=10 0 0;0 -10 0;0 0 10;L=0 0 0;-2 0 0;-1 -3 0;U=0 -3 -1;0 0 -3; 0 0 0;w=1.1;for i=1:9 x_sor(:,i)=inv(D./w-L)*(1/w-1)*D+U)*x0+inv(D./w-L)*b; x0=x_sor(:,i);endx_sor P131 4-7A=3 2;2 6;b=2 -8;x0=-2 -2; for i=1:8 r0=b-A*x0; a0=r0*r0/(r0*A*r0); x1=x0+a0*r0; x(:,i)=x1; x0=x1;endxP157, 5-6a=3,1;1,3;v0=0,1;for i=1:10v=a*v0;max(v)V=v/max(v)Eig=(V*a*V)/(V*V)v0=V;endP160, 5-8a=3,1;1,3;v0=0,1;for i=1:10v=inv(a)*v0;max(v)V=v/max(v)v0=V;end%P202 example6-6A=1 1 1 1 1;1.00 1.25 1.5 1.75 2.00;f=1.6292 1.7561 1.8764 2.0082 2.1353; x=inv(A*A)*A*f %P205 example 6-8Sigma1=sqrt(sum(A(:,1).2);v1=A(:,1)+Sigma1*1 0 0 0 0;A1=A(:,1)-2*(v1*A(:,1)./(v1*v1)*v1A2=A(:,2)-2*(v1*A(:,2)./(v1*v1)*v1f2=f-2*(v1*f)./(v1*v1)*v1AA=A1,A2 A2(1)=0;Sigma2=-sqrt(sum(A2.2)v2=A2+Sigma2*0 1 0 0 0;f3=f2-2*(v2*f2)./(v2*v2)*v2A3=A2-2*(v2*A2)./(v2*v2)*v2A3(1)=AA(1,2)AAA=A1,A3 x=inv(AAA(1:2,1:2)*f3(1:2,1)P218, 6-11 x=0.4 0.55 0.65 0.8 0.9 1.05;f=0.41075 0.57815 0.69675 0.88811 1.02652 1.25382; diff_1=(f(2:6)-f(1:5)./(x(2:6)-x(1:5) diff_2=(diff_1(2:5)-diff_1(1:4)./(x(3:6)-x(1:4) diff_3=(diff_2(2:4)-diff_2(1:3)./(x(4:6)-x(1:3) diff_4=(diff_3(2:3)-diff_3(1:2)./(x(5:6)-x(1:2) diff_5=(diff_4(2)-diff_4(1)./(x(6)-x(1)P230, 6-13 x=1:6;y=16 18 21 17 15 12;xi=0.75:0.05:6.25;p=polyfit(x,y,5);v1=polyval(p,xi);v2=interp1(x,y,xi,linear);v3=interp1(x,y,xi,pchip);v4=interp1(x,y,xi,spline);subplot(2,2,1);plot(x,y,o,xi,v1,-);subplot(2,2,2);plot(x,y,o,xi,v2,-);subplot(2,2,3);plot(x,y,o,xi,v3,-);subplot(2,2,4);plot(x,y,o,xi,v4,-);%P249 example 7-4x=0:1/8:1;f=sin(x)./x;f(1)=1; % T n=8;h=(1-0)/n;T=1/2*h*(f(1)+f(9)+2*sum(f(2:8) %Sn=4;h=(1-0)/n;S=h/6*(f(1)+f(9)+4*sum(f(2:2:8)+2*sum(f(3:2:7) %P255 example 7-5%第二列T01=1/2*(f(1)+f(9)T02=T01/2+1/2*(f(5)T03=T02/2+1/4*(f(3)+f(7) %第三列T11=1/3*(4*T02-T01)T12=1/3*(4*T03-T02) %第四列T21=1/15*(16*T12-T11)%P255 example 7-6format longT0=0.5;h=1;for i=1:5 a=1:2:2i-1; x=a/(2i); sum1=sum(x.(3/2); T=1/2*T0+h/2*sum1 T0=T; TT(i)=T; h=h/2; endTT=0.5,TTTT1=1/3*(4*TT(2:6)-TT(1:5)TT2=1/15*(16*TT1(2:5)-TT1(1:4)TT3=1/(43-1)*(43*TT2(2:4)-TT2(1:3)TT4=1/(44-1)*(44*TT3(2:3)-TT3(1:2)TT5=1/(45-1)*(45*TT4(2)-TT4(1)P267, 7-10for i=1:12 h=10(-i); x=2; Dc=(sqrt(x+h)-sqrt(x-h)/(2*h)endP270, 7-11 h=0.1;x=0.5;x1=x+h;x2=x-h;Dc1=(x12*exp(-x1)-x22*exp(-x2)/(2*h); h=0.05;x=0.5;x1=x+h;x2=x-h;Dc2=(x12*exp(-x1)-x22*exp(-x2)/(2*h); h=0.025;x=0.5;x1=x+h;x2=x-h;Dc3=(x12*exp(-x1)-x22*exp(-x2)/(2*h); D1c1=1/3*(4*Dc2-Dc1)D1c2=1/3*(4*Dc3-Dc2) D2c1=1/15*(16*D1c2-D1c1)P288, 8-7format longy1=0.4;h=0.1;