end p=S(1);
利用上述程序计算p(x)在x?2邻域附近的值具体见下表：
x p(x) 1.95 -9.6634e-13 1.96 -1.1255e-11 1.97 2.8422e-12 1.98 7.3896e-12 1.99 5.3433e-12 2.00 0 x p(x) 2.01 -3.7517e-12 2.02 -7.7875e-12 2.03 4.4338e-12 2.04 -3.6380e-12 2.05 5.1159e-12 当x?[1.95,20.5]时，p(x)的图像如下： 画图程序如下： function huatu(A,x) [~,n]=size(x); for i=1:n
y(i)=qinjiushao(A,x(i)); end plot(x,y); 程序运行如下： >> x=1.95:0.01:20.5;
>> A=[-512 2304 -4608 5376 -4032 2016 -672 144 -18 1]; >> huatu(A,x)
3x 10112.521.5p(x)10.50-0.50510x152025
4：编制计算机给定矩阵A的LU分解和PLU分解的通用程序，然后利用你编写的程序完成下面两个计算任务：
考虑
?1??1?A??????1???10?1?????01??Rn?n
?11??11??0??????1??1自己取定x?Rn，并计算b?Ax。然后用你编制的不选主元的Gauss消去法求解
?。对n从5到30估计计算解的精度。 该方程组，记你计算出的解为x（2）对n从5到30计算出其逆矩阵。 LU分解的通用程序： function LU(A) [m,n]=size(A); L=zeros(m,n); U=zeros(m,n); for j=1:n
U(1,j)=A(1,j); L(j,j)=1;
end for j=2:m
L(j,1)=A(j,1)/U(1,1); end for i=2:n for j=i:n sum=0; for k=1:i-1
sum=sum+L(i,k)*U(k,j); end
U(i,j)=A(i,j)-sum; end for j=i+1:n sum=0; for k=1:i-1
sum=sum+L(j,k)*U(k,i); end
L(j,i)=(A(j,i)-sum)/U(i,i); end end disp('L ='); disp(L); disp('U ='); disp(U);
PLU分解的通用程序： function PLU(A)
[~,n]=size(A); Ip=1:n; for k=1:n-1
[~,r]=max(abs(A(k:n,k))); r=r+(k-1); if r>k
A([k,r],:)=A([r,k],:); Ip([k,r])=Ip([r,k]); end for p=k+1:n mu=A(p,k)/A(k,k); A(p,k)=mu;
A(p,k+1:n)=A(p,k+1:n)-mu*A(k,k+1:n); end end p=eye(n,n); P=zeros(n,n); for i=1:n
P(i,1:n)=p(Ip(i),1:n); end
L=tril(A,-1)+eye(n) ; U=triu(A); disp('P=') disp(P); disp('L=') disp(L);
disp('U='); disp(U);
我选取b?[123?n]Tn?[5,30],n?N*，则我们可以计算出其精确解
2n?1?12n?2?1?1x?[?n?1?n?2?2222n?1T]2n?1n?[5,30]n?N*
实现求解的程序如下： function INV(n) A=ones(n); for i=2:n for j=1:i-1
A(i,j)=-1*A(i,j); end for j=i:n-1 A(i-1,j)=0; end end for i=1:n b(i)=i;
X(i)=-(2^(n-i)-1)/(2^(n-i)); end
X(n)=(2^n-1)/(2^(n-1)); [L,U]=LU(A); x1=inv(U)*inv(L)*b' disp(x1);
我们利用上述的程序计算出的结果如下表：
n x1,x2,x3,?,xn