清华大学数值分析A第一次作业
数值分析作业及参考答案
数值分析第一次作业及参考答案1. 设212S gt =，假定g 是准确的，而对t 的测量有0.1±秒的误差，证明当t 增加时S 的绝对误差增加，而相对误差却减少。
解：2**22211()0.122()0.10.2()1122,(),().r r e S S S gt gt gt e S gt e S t gt gt t e S e S =-=-====∴↑↑↓2. 设2()[,]f x C a b ∈且()()0f a f b ==，求证2''1max ()()max ().8a x ba xb f x b a f x ≤≤≤≤≤-解：由112,0),(,0)()()0()00.a b L x l x l x =⨯+⨯=（两点线性插值 插值余项为"111()()()()()()[,]2R x f x L x f x a x b a b ξξ=-=--∈ [,].x a b ∴∀∈有12211()()"()()()max "()[()()]221()()1max "()[]()max "().228a x ba xb a x b f x R x f x a x b f x x a b x x a b x f x b a f x ξ≤≤≤≤≤≤==--≤---+-≤=-21max ()()max "()8a xb a x b f x b a f x ≤≤≤≤∴≤-3. 已测得函数()y f x =的三对数据：（0，1），（－1，5），（2，－1），（1）用Lagrange 插值求二次插值多项式。
（2）构造差商表。
（3）用Newton 插值求二次插值多项式。
解：(1)Lagrange 插值基函数为0(1)(2)1()(1)(2)(01)(02)2x x l x x x +-==-+-+-同理 1211()(2),()(1)36l x x x l x x x =-=+ 故2202151()()(1)(2)(2)(1)23631i i i p x y l x x x x x x x x x =-==-+-+-++=-+∑（2）令0120,1,2x x x ==-=，则一阶差商、二阶差商为0112155(1)[,]4,[,]20(1)12f x x f x x ---==-==-----0124(2)[,,]102f x x x ---==-22()1(4)(0)1*(0)(1)31P x x x x x x =+--+-+=-+4. 在44x -≤≤上给出()xf x e =的等距节点函数表，若用二次插值求x e 的近似值，要使截断误差不超过610-，问使用函数表的步长h 应取多少？解：()40000(),(),[4,4],,,, 1.x k x f x e f x e e x x h x x h x x th t ==≤∈--+=+≤考察点及(3)200044343()()[(()]()[()]3!(1)(1)(1)(1)3!3!.(4,4).6f R x x x h x x x x h t t t e t h th t h e h e ξξ=----+-+≤+⋅⋅-=≤∈-则436((1)(1)100.006.t t t h --+±<< 在点 得5. 求2()f x x =在［a,b ］上的分段线性插值函数()h I x ，并估计误差。
数值分析大作业一
数值分析大作业一一、算法设计方案1、求λ1和λ501的值：思路：采用幂法求出按模最大特征值λmax，该值必为λ1或λ501，若λmax小于0，则λmax=λ1；否则λmax=λ501。
再经过原点平移，使用幂法迭代出矩阵A-λmax I的特征值，此时求出的按模最大特征值即为λ1和λ501的另一个值。
2、求λs的值：采用反幂法求出按模最小的特征值λmin即为λs，其中的方程组采用LU分解法进行求解。
3、求与μk最接近的特征值：对矩阵A采用带原点平移的反幂法求解最小特征值，其中平移量为：μk。
4、A的条件数cond(A)=| λmax/λmin|；5、A的行列式的值：先将A进行LU分解，再求U矩阵对角元素的乘积即为A 行列式的值。
二、源程序#include<iostream>#include<iomanip>#include<math.h>#define N 501#define E 1.0e-12 //定义精度常量#define r 2#define s 2using namespace std;double a[N];double cc[5][N];void init();double mifa();double fmifa();int max(int aa,int bb);int min(int aa,int bb);int max_3(int aa,int bb,int cc);void LU();void main(){double a1,a2,d1,d501=0,ds,det=1,miu[39],lamta,cond;int i,k;init();/*************求λ1和λ501********************/a1=mifa();if(a1<0)d1=a1; //若小于0则表示λ1的值elsed501=a1; //若大于0则表示λ501的值for(i=0;i<N;i++)a[i]=a[i]-a1;a2=mifa()+a1;if(a2<0)d1=a2; //若小于0则表示λ1的值elsed501=a2; //若大于0则表示λ501的值cout<<"λ1="<<setiosflags(ios::scientific)<<setprecision(12)<<d1<<"\t";cout<<"λ501="<<setiosflags(ios::scientific)<<setprecision(12)<<d501<<endl;/**************求λs*****************/init();ds=fmifa();cout<<"λs="<<setiosflags(ios::scientific)<<setprecision(12)<<ds<<endl;/**************求与μk最接近的特征值λik**************/cout<<"与μk最接近的特征值λik:"<<endl;for(k=0;k<39;k++){miu[k]=d1+(k+1)*(d501-d1)/40;init();for(i=0;i<N;i++)a[i]=a[i]-miu[k];lamta=fmifa()+miu[k];cout<<"λi"<<k+1<<"\t\t"<<setiosflags(ios::scientific)<<setprecision(12)<<lamta<<en dl;}/**************求A的条件数**************/cout<<"矩阵A的条件式";cond=abs(max(abs(d1),abs(d501))/ds);cout<<"cond="<<setiosflags(ios::scientific)<<setprecision(12)<<cond<<endl;/**************求A的行列式**************/cout<<"矩阵A的行列式";init();LU();for(i=0;i<N;i++){det*=cc[2][i];}cout<<"det="<<setiosflags(ios::scientific)<<setprecision(12)<<det<<endl;system("pause");}/**************初始化函数，给a[N]赋值*************/void init(){int i;for(i=1;i<=501;i++)a[i-1]=(1.64-0.024*i)*sin((double)(0.2*i))-0.64*exp((double)(0.1/i)); }/**************幂法求最大绝对特征值**************/double mifa(){int i,k=0;double u[N],y[N]={0},b=0.16,c=-0.064,Beta_=0,error;for(i=0;i<501;i++)u[i]=1; //令u[N]=1for(k=1;k<2000;k++) //控制最大迭代次数为2000{/***求y(k-1)***/double sum_u=0,gh_sum_u;for(i=0;i<N;i++){sum_u+=u[i]*u[i]; }gh_sum_u=sqrt(sum_u);for(i=0;i<N;i++){y[i]=u[i]/gh_sum_u;}/****求新的uk****/u[0]=a[0]*y[0]+b*y[1]+c*y[2];u[1]=b*y[0]+a[1]*y[1]+b*y[2]+c*y[3]; //前两列和最后两列单独拿出来求中D间的循环求for(i=2;i<N-2;i++){u[i]=c*y[i-2]+b*y[i-1]+a[i]*y[i]+b*y[i+1]+c*y[i+2];}u[N-2]=c*y[N-4]+b*y[N-3]+a[N-2]*y[N-2]+b*y[N-1];u[N-1]=c*y[N-3]+b*y[N-2]+a[N-1]*y[N-1];/***求beta***/double Beta=0;for(i=0;i<N;i++){Beta+=y[i]*u[i];}//cout<<"Beta"<<k<<"="<<Beta<<"\t"; 输出每次迭代的beta /***求误差***/error=abs(Beta-Beta_)/abs(Beta);if(error<=E) //若迭代误差在精度水平内则可以停止迭代{return Beta;} //控制显示位数Beta_=Beta; //第个eta的值都要保存下来，为了与后个值进行误差计算 }if(k==2000){cout<<"error"<<endl;return 0;} //若在最大迭代次数范围内都不能满足精度要求说明不收敛}/**************反幂法求最小绝对特¬征值**************/double fmifa(){int i,k,t;double u[N],y[N]={0},yy[N]={0},b=0.16,c=-0.064,Beta_=0,error;for(i=0;i<501;i++)u[i]=1; //令u[N]=1for(k=1;k<2000;k++){double sum_u=0,gh_sum_u;for(i=0;i<N;i++){sum_u+=u[i]*u[i]; }gh_sum_u=sqrt(sum_u);for(i=0;i<N;i++){y[i]=u[i]/gh_sum_u;yy[i]=y[i]; //用重新赋值，避免求解方程组的时候改变y的值}/****LU分解法解方程组Au=y，求新的***/LU();for(i=2;i<=N;i++){double temp_b=0;for(t=max(1,i-r);t<=i-1;t++)temp_b+=cc[i-t+s][t-1]*yy[t-1];yy[i-1]=yy[i-1]-temp_b;}u[N-1]=yy[N-1]/cc[s][N-1];for(i=N-1;i>=1;i--){double temp_u=0;for(t=i+1;t<=min(i+s,N);t++)temp_u+=cc[i-t+s][t-1]*u[t-1];u[i-1]=(yy[i-1]-temp_u)/cc[s][i-1];}double Beta=0;for(i=0;i<N;i++){Beta+=y[i]*u[i];}error=abs(Beta-Beta_)/abs(Beta);if(error<=E){return (1/Beta);}Beta_=Beta;}if(k==2000){cout<<"error"<<endl;return 0;} }/**************求两数最大值的子程序**************/int max(int aa,int bb){return(aa>bb?aa:bb);}/**************求两数最小值的子程序**************/int min(int aa,int bb){return(aa<bb?aa:bb);}/**************求三数最大值的子程序**************/int max_3(int aa,int bb,int cc){ int tt;if(aa>bb)tt=aa;else tt=bb;if(tt<cc) tt=cc;return(tt);}/**************LU分解**************/void LU(){int i,j,k,t;double b=0.16,c=-0.064;/**赋值压缩后矩阵cc[5][501]**/for(i=2;i<N;i++)cc[0][i]=c;for(i=1;i<N;i++)cc[1][i]=b;for(i=0;i<N;i++)cc[2][i]=a[i];for(i=0;i<N-1;i++)cc[3][i]=b;for(i=0;i<N-2;i++)cc[4][i]=c;for(k=1;k<=N;k++){for(j=k;j<=min(k+s,N);j++){double temp=0;for(t=max_3(1,k-r,j-s);t<=k-1;t++)temp+=cc[k-t+s][t-1]*cc[t-j+s][j-1];cc[k-j+s][j-1]=cc[k-j+s][j-1]-temp;}//if(k<500){for(i=k+1;i<=min(k+r,N);i++){double temp2=0;for(t=max_3(1,i-r,k-s);t<=k-1;t++)temp2+=cc[i-t+s][t-1]*cc[t-k+s][k-1];cc[i-k+s][k-1]=(cc[i-k+s][k-1]-temp2)/cc[s][k-1];}}}}三、程序结果。
(完整版)数值分析第一次作业
问题1：20.给定数据如下表：试求三次样条插值S(x)，并满足条件 (1)S`(0.25)=1.0000，S`(0.53)=0.6868； （2）S ’’(0.25)=S ’’(0.53)=0。
分析：本问题是已知五个点，由这五个点求一三次样条插值函数。
边界条件有两种，（1）是已知一阶倒数，（2）是已知自然边界条件。
对于第一种边界(已知边界的一阶倒数值)，可写出下面的矩阵方程。
⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎣⎡=⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎣⎡⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎣⎡432104321034322110d M M M M M 200020000020022d d d d λμμλμλμλ其中μj =j1-j 1-j h h h +，λi=j1-j j h h h +，dj=6f[x j-1,x j ,x j+1]， μn =1，λ0=1对于第一种边界条件d 0=0h 6（f[x 0,x 1]-f 0`），d n =1-n h 6（f`n-f `[x n-1,x n ]） 解：由matlab 计算得：由此得矩阵形式的线性方程组为：⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎣⎡=⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎣⎡⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎣⎡ 2.1150-2.4286-3.2667-4.3143-5.5200-M M M M M 25714.00001204286.000004000.026000.0006429.023571.0001243210解得 M 0=-2.0286;M 1=-1.4627;M 2= -1.0333; M 3= -0.8058; M 4=-0.6546S(x)=⎪⎪⎩⎪⎪⎨⎧∈-+-+-∈-+-+-∈-+-+-∈-+-+-]53.0,45.0[x 5.40x 9.1087x 35.03956.8.450-x 1.3637-x .5301.67881- ]45.0,39.0[x 9.30x 11.188x 54.010.418793.0-x 2.2384-x .450(2.87040-]39.0,30.0[x 03.0x 6.9544x 9.30 6.107503.0-x 1.9136-x .3902.708779-]30.0,25.0[x 5.20x 10.9662x 0.3010.01695.20-x 4.8758-x .3006.76209-33333333），（）（）（）（），（）（）（）），（）（）（）（），（）（）（）（Matlab 程序代码如下：function tgsanci(n,s,t) %n代表元素数，s,t代表端点的一阶导。
清华大学高等数值分析作业李津1——矩阵基础
20130917题目求证：在矩阵的LU 分解中，111n n Tn ij i j j i j L I e e α-==+⎛⎫=- ⎪⎝⎭∑∑证明：在高斯消去过程中，假设0jj a ≠ ，若a=0，可以通过列变换使得前面的条件成立，这里不考虑这种情况。
对矩阵A 进行LU 分解，()()()()()1111111L M n M M M n ---=-=••-………… ，其中()1n Tn ij i j i j M j I e e α=+⎛⎫=+ ⎪⎝⎭∑ ，i e 、j e 为n 维线性空间的自然基。
()M j 是通过对单位阵进行初等变换得到，通过逆向的变换则可以得到单位阵，由此很容易得到()M j 的逆矩阵为1n T n ij i j i j I e e α=+⎛⎫- ⎪⎝⎭∑。
故111n n T n ij i j n j i j L I e e I α-==+⎛⎫⎛⎫=- ⎪ ⎪ ⎪⎝⎭⎝⎭∏∑上式中的每一项均是初等变换，从右向左乘，则每乘一次相当于对右边的矩阵进行一次向下乘法叠加的初等变换。
由于最初的矩阵为单位阵，变换从右向左展开，因而每一次变换不改变已经更新的数据，既该变换是从右向左一列一列更新数据，故11nn Tn ij i j j i j L I e e α==+⎛⎫=- ⎪⎝⎭∑∑。
数学证明：1n Tij i j i j e e α=+⎛⎫ ⎪⎝⎭∑具有,000n j j A -⎛⎫ ⎪⎝⎭ 和1,1000n j n j B -+-+⎛⎫⎪⎝⎭ 的形式，且有+1,-11,10000=000n j j n j n j A B --+-+⎛⎫⎛⎫⎪⎪⎝⎭⎝⎭而11n n T ij i j j k i j e e α-==+⎛⎫ ⎪⎝⎭∑∑具有1,1000n k n k B -+-+⎛⎫⎪⎝⎭的形式，因此： 1311111211121==n n n n n n T T T n ij i j n ij i j n ik i k j i j j i j k n i k n n T n i i n ik i i i k L I e e I e e I e e I e e I e ααααα---==+==+=-=+==+⎡⎤⎛⎫⎛⎫⎛⎫⎛⎫⎛⎫⎛⎫=---⎢⎥ ⎪ ⎪ ⎪ ⎪ ⎪⎪ ⎪ ⎪⎝⎭⎝⎭⎢⎥⎝⎭⎝⎭⎝⎭⎝⎭⎣⎦⎛⎫⎛⎫⎛⎫=-- ⎪ ⎪ ⎝⎭⎝⎝⎭∏∑∏∑∑∑∑∑……11211n n n T Tk n ik i kk k i k e I e e α--===+⎛⎫⎛⎫=- ⎪⎪ ⎪⎭⎝⎭⎝⎭∑∑∑#20130924题目一问：能否用逐次householder 相似变换变实矩阵A 为上三角矩阵，为什么？解：不能用逐次householder 相似变换变A 为上三角矩阵，原因如下：A 记作：()12=,,n A a a a ……， ，存在householder 阵1H s.t. 1111H a e α= ，则()()()111111111111111111111,,,0T Th H AH H a A H e H A H e H A H h H A H ααα⎛⎫'''=== ⎪⎪'⎝⎭⎛⎫''=+ ⎪ ⎪⎝⎭11H A H ''第一列的元素不能保证为1e 的倍数，故无法通过householder 变换实现上三角化。
清华大学高等数值分析 第一次实验作业
10
-10
0
100
200
300
400
500
600
700
800
900
迭代次数
图9
m=100时，Lanczos法求解Ax=b的收敛曲线
高等数值分析实验作业一
10
4
Lanzcos 算法的收敛曲线 (阶数 n=1002）
10
2
10
0
||rk||/||b||
10
-2
10
-4
10
-6
10
-8
10
-10
0
200
迭代次数
图12 m=10时，Minres法求解Ax=b的收敛曲线
10
2
Minres 算法的收敛曲线 (阶数 n=1002）
10
0
10
-2
||rk||/||b||
10
-4
10
-6
10
-8
10
-10
0
100
200
300
400
500
600
700
迭代次数
图13
10
2
m=50时，Minres法求解Ax=b的收敛曲线
10
0
Lanzcos 算法的收敛曲线 (阶数 n=1002)
m=10 m=50 m=100 m=400 m=800
10
-2
10
-4
||rk||/||b||
10
-6
10
-8
10
-10
10
-12
0
2
4
6
8
10
12
14
16
清华大学高等数值分析实验设计及答案
高等数值分析实验一工物研13 成彬彬2004310559一．用CG，Lanczos和MINRES方法求解大型稀疏对称正定矩阵Ax=b作实验中，A是利用A= sprandsym(S,[],rc,3)随机生成的一个对称正定阵，S是1043阶的一个稀疏阵A= sprandsym(S,[],0.01,3);检验所生成的矩阵A的特征如下：rank(A-A')=0 %即A=A’，A是对称的；rank(A)=1043 %A满秩cond(A)= 28.5908 %A是一个“好”阵1．CG方法利用CG方法解上面的线性方程组[x,flag,relres,iter,resvec] = pcg(A,b,1e-6,1043);结果如下：Iter=35，表示在35步时已经收敛到接近真实xrelres= norm(b-A*x)/norm(b)= 5.8907e-007为最终相对残差绘出A的特征值分布图和收敛曲线：S=svd(A); %绘制特征值分布subplot(211)plot(S);title('Distribution of A''s singular values');;xlabel('n')ylabel('singular values')subplot(212); %绘制收敛曲线semilogy(0:iter,resvec/norm(b),'-o');title('Convergence curve');xlabel('iteration number');ylabel('relative residual');得到如下图象：为了观察CG方法的收敛速度和A的特征值分布的关系，需要改变A的特征值：(1).研究A的最大最小特征值的变化对收敛速度的影响在A的构造过程中，通过改变A= sprandsym(S,[],rc,3)中的参数rc(1/rc为A的条件数)，可以达到改变A的特征值分布的目的：通过改变rc=0.1，0.0001得到如下两幅图以上三种情况下，由收敛定理2.2.2计算得到的至多叠代次数分别为：48，14和486，由于上实验结果可以看出实际叠代次数都比上限值要小较多。
(完整版)数值分析第一次作业
问题1：20.给定数据如下表：试求三次样条插值S(x)，并满足条件 (1)S`(0.25)=1.0000，S`(0.53)=0.6868； （2）S ’’(0.25)=S ’’(0.53)=0。
分析：本问题是已知五个点，由这五个点求一三次样条插值函数。
边界条件有两种，（1）是已知一阶倒数，（2）是已知自然边界条件。
对于第一种边界(已知边界的一阶倒数值)，可写出下面的矩阵方程。
⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎣⎡=⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎣⎡⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎣⎡432104321034322110d M M M M M 200020000020022d d d d λμμλμλμλ其中μj =j1-j 1-j h h h +，λi=j1-j j h h h +，dj=6f[x j-1,x j ,x j+1]， μn =1，λ0=1对于第一种边界条件d 0=0h 6（f[x 0,x 1]-f 0`），d n =1-n h 6（f`n-f `[x n-1,x n ]） 解：由matlab 计算得：由此得矩阵形式的线性方程组为：⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎣⎡=⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎣⎡⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎣⎡ 2.1150-2.4286-3.2667-4.3143-5.5200-M M M M M 25714.00001204286.000004000.026000.0006429.023571.0001243210解得 M 0=-2.0286;M 1=-1.4627;M 2= -1.0333; M 3= -0.8058; M 4=-0.6546S(x)=⎪⎪⎩⎪⎪⎨⎧∈-+-+-∈-+-+-∈-+-+-∈-+-+-]53.0,45.0[x 5.40x 9.1087x 35.03956.8.450-x 1.3637-x .5301.67881- ]45.0,39.0[x 9.30x 11.188x 54.010.418793.0-x 2.2384-x .450(2.87040-]39.0,30.0[x 03.0x 6.9544x 9.30 6.107503.0-x 1.9136-x .3902.708779-]30.0,25.0[x 5.20x 10.9662x 0.3010.01695.20-x 4.8758-x .3006.76209-33333333），（）（）（）（），（）（）（）），（）（）（）（），（）（）（）（Matlab 程序代码如下：function tgsanci(n,s,t) %n代表元素数，s,t代表端点的一阶导。
数值分析第一次大作业
《数值分析》计算实习报告第一题院系：机械工程及自动化学院_学号： _____姓名： _ ______2017年11月7日一、算法设计方案1、求λ1，λ501和λs 的值1)利用幂法计算出矩阵A 按模最大的特征值，设其为λm 。
2)令矩阵B =A −λm I （I 为单位矩阵），同样利用幂法计算出矩阵B 按模最大的特征值λm ′。
3)令λm ′′=λm ′+λm 。
由计算过程可知λm 和λm ′′分别为矩阵A 所有特征值按大小排序后，序列两端的值。
即，λ1=min{λm ,λm ′′}，λ501=max{λm ,λm ′′}。
4) 利用反幂法计算λs 。
其中，反幂法每迭代一次都要求解线性方程组1k k Au y -=，由于矩阵A 为带状矩阵，故可用三角分解法解带状线性方程组的方法求解得到k u 。
2、求A 的与数μk =λ1+k λ501−λ140最接近的特征值λi k (k =1,2, (39)1) 令矩阵D k =A −μk I ，利用反幂法计算出矩阵D k 按模最小的特征值λi k ′,则λi k =λi k ′+μk 。
3、求A 的（谱范数）条件数cond(A )2和行列式det A1) cond(A)2=|λm λs |，前文已算出m λ和s λ，直接带入即可。
2) 反幂法计算λs 时，已经对矩阵A 进行过Doolittle 分解，得到A=LU 。
而L 为对角线上元素全为1的下三角矩阵，U 为上三角矩阵，可知det 1L =，5011det ii i U u ==∏，即有5011det det det ii i A L U u ====∏。
最后，为节省存储量，需对矩阵A 进行压缩，将A 中带内元素存储为数组C [5][501]。
二、源程序代码#include<windows.h>#include<iostream>#include<iomanip>#include<math.h>using namespace std;#define N 501#define K 39#define r 2#define s 2#define EPSI 1.0e-12//求两个整数中的最大值int int_max2(int a, int b){return(a>b ? a : b);}//求两个整数中的最小值int int_min2(int a, int b){return(a<b ? a : b);}//求三个整数中的最大值int int_max3(int a, int b, int c){int t;if (a>b)t = a;else t = b;if (t<c) t = c;return(t);}//定义向量内积double dianji(double x[], double y[]) {double sum = 0;for (int i = 0; i<N; i++)sum = sum + x[i] * y[i];return(sum);}//计算两个数之间的相对误差double erro(double lamd0, double lamd1){double e, d, l;e = fabs(lamd1 - lamd0);d = fabs(lamd1);l = e / d;return(l);}//矩阵A的压缩存储初始化成Cvoid init_c(double c[][N]){int i, j;for (i = 0; i<r + s + 1; i++)for (j = 0; j<N; j++)if (i == 0 || i == 4)c[i][j] = -0.064;else if (i == 1 || i == 3)c[i][j] = 0.16;elsec[i][j] = (1.64 - 0.024*(j + 1))*sin(0.2*(j + 1)) - 0.64*exp(0.1 / (j + 1)); }//矩阵复制void fuzhi_c(double c_const[][N], double c[][N]){int i, j;for (i = 0; i<r + s + 1; i++)for (j = 0; j<N; j++)c[i][j] = c_const[i][j];}//LU三角分解void LUDet_c(double c_const[][N], double c_LU[][N]){double sum;int k, i, j;fuzhi_c(c_const, c_LU);for (k = 1; k <= N; k++){for (j = k; j <= int_min2(k + s, N); j++){sum = 0;for (i = int_max3(1, k - r, j - s); i <= k - 1; i++)sum += c_LU[k - i + s][i - 1] * c_LU[i - j + s][j - 1];c_LU[k - j + s][j - 1] -= sum;}for (j = k + 1; j <= int_min2(k + r, N); j++){sum = 0;for (i = int_max3(1, j - r, k - s); i <= k - 1; i++)sum += c_LU[j - i + s][i - 1] * c_LU[i - k + s][k - 1];c_LU[j - k + s][k - 1] = (c_LU[j - k + s][k - 1] - sum) / c_LU[s][k - 1];}}}//三角分解法解带状线性方程组void jiefc(double c_const[][N], double b_const[], double x[]){int i, j;double b[N], c_LU[r + s + 1][N], sum;for (i = 0; i<N; i++)b[i] = b_const[i];LUDet_c(c_const, c_LU);for (i = 2; i <= N; i++){sum = 0;for (j = int_max2(i - 2, 1); j <= i - 1; j++)sum += c_LU[i - j + 2][j - 1] * b[j - 1];b[i - 1] -= sum;}x[N - 1] = b[N - 1] / c_LU[2][N - 1];for (i = N - 1; i >= 1; i--){sum = 0;for (j = i + 1; j <= int_min2(i + 2, N); j++)sum += c_LU[i - j + 2][j - 1] * x[j - 1];x[i - 1] = (b[i - 1] - sum) / c_LU[2][i - 1];}}//幂法求按模最大特征值double mifa_c(double c_const[][N]){double u[N], y[N];double sum, length_u, beta0, beta1;int i, j;for (i = 0; i<N; i++)//迭代初始向量u[i] = 0.5;length_u = sqrt(dianji(u, u));for (i = 0; i<N; i++)y[i] = u[i] / length_u;for (i = 1; i <= N; i++){sum = 0;for (j = int_max2(i - 2, 1); j <= int_min2(i + 2, N); j++)sum = sum + c_const[i - j + 2][j - 1] * y[j - 1];u[i - 1] = sum;}beta1 = dianji(u, y);do{beta0 = beta1;length_u = sqrt(dianji(u, u));for (i = 0; i<N; i++)y[i] = u[i] / length_u;for (i = 1; i <= N; i++){sum = 0;for (j = int_max2(i - 2, 1); j <= int_min2(i + 2, N); j++)sum = sum + c_const[i - j + 2][j - 1] * y[j - 1];u[i - 1] = sum;}beta1 = dianji(u, y);} while (erro(beta0, beta1) >= EPSI);return(beta1);}//反幂法求按模最小特征值double fmifa_c(double c_const[][N]){double u[N], y[N];double length_u, beta0, beta1;int i;for (i = 0; i<N; i++)//迭代初始向量u[i] = 0.5;length_u = sqrt(dianji(u, u));for (i = 0; i<N; i++)y[i] = u[i] / length_u;jiefc(c_const, y, u);beta1 = dianji(y, u);do{beta0 = beta1;length_u = sqrt(dianji(u, u));for (i = 0; i<N; i++)y[i] = u[i] / length_u;jiefc(c_const, y, u);beta1 = dianji(y, u);} while (erro(beta0, beta1) >= EPSI);beta1 = 1 / beta1;return(beta1);}//计算lamd_1、lamd_501、lamd_svoid calculate1(double c_const[][N], double &lamd_1, double &lamd_501, double &lamd_s) {int i;double lamd_mifa0, lamd_mifa1, c[r + s + 1][N];lamd_mifa0 = mifa_c(c_const);fuzhi_c(c_const, c);for (i = 0; i<N; i++)c[2][i] = c[2][i] - lamd_mifa0;lamd_mifa1 = mifa_c(c) + lamd_mifa0;if (lamd_mifa0<lamd_mifa1){lamd_1 = lamd_mifa0;lamd_501 = lamd_mifa1;}else{lamd_501 = lamd_mifa0;lamd_1 = lamd_mifa1;}lamd_s = fmifa_c(c_const);}//平移+反幂法求最接近u_k的特征值void calculate2(double c_const[][N], double lamd_1, double lamd_501, double lamd_k[]){int i, k;double c[r + s + 1][N], h, temp;temp = (lamd_501 - lamd_1) / 40;for (k = 1; k <= K; k++){h = lamd_1 + k*temp;fuzhi_c(c_const, c);for (i = 0; i<N; i++)c[2][i] = c[2][i] - h;lamd_k[k - 1] = fmifa_c(c) + h;}}//计算cond(A)和det(A)void calculate3(double c_const[][N], double lamd_1, double lamd_501, double lamd_s, double &cond_A, double &det_A){int i;double c_LU[r + s + 1][N];if (fabs(lamd_1)>fabs(lamd_501))cond_A = fabs(lamd_1 / lamd_s);elsecond_A = fabs(lamd_501 / lamd_s);LUDet_c(c_const, c_LU);det_A = 1;for (i = 0; i<N; i++)det_A *= c_LU[2][i];}//*主程序*//int main(){int i, count = 0;double c_const[5][N], lamd_k[K];double lamd_1, lamd_501, lamd_s;double cond_A, det_A;//设置白背景黑字system("Color f0");//矩阵A压缩存储到c[5][501]init_c(c_const);cout << setiosflags(ios::scientific) << setiosflags(ios::right) << setprecision(12) << endl;//计算lamd_1、lamd_501、lamd_scalculate1(c_const, lamd_1, lamd_501, lamd_s);cout << " 矩阵A的最小特征值：λ1 = " << setw(20) << lamd_1 << endl;cout << " 矩阵A的最大特征值：λ501 = " << setw(20) << lamd_501 << endl;cout << " 矩阵A的按模最小的特征值：λs = " << setw(20) << lamd_s << endl;//求最接近u_k的特征值calculate2(c_const, lamd_1, lamd_501, lamd_k);cout << endl << " 与数u_k最接近的特征值：" << endl;for (i = 0; i<K; i++){cout << " λ_ik_" << setw(2) << i + 1 << " = " << setw(20) << lamd_k[i] << " ";count++;if (count == 2){cout << endl;count = 0;}}//计算cond_A和det_Acalculate3(c_const, lamd_1, lamd_501, lamd_s, cond_A, det_A);cout << endl << endl;cout << " 矩阵A的条件数：cond(A) = " << setw(20) << cond_A << endl;cout << " 矩阵A的行列式的值：det(A) = " << setw(20) << det_A << endl << endl;return 0;}三，计算结果四，分析初始向量选择对计算结果的影响当选取初始向量0(1,1,,1)Tu=时，计算的结果如下：此结果即为上文中的正确计算结果。
数值分析第五版习题答案清华大学出版社
若，则
则
则法方程组为
从而解得
故
均方误差为
21。在某佛堂反应中，由实验得分解物浓度与时间关系如下：
时间
0 5 10 15 20 25 30 35 40 45 50 55
浓度
0 1.27 2.16 2.86 3.44 3.87 4.15 4.37 4.51 4.58 4.62 4.64
用最小二乘法求。
计算到。若取（5位有效数字），试问计算将有多大误差？
解：
……
依次代入后，有
即，
若取,
的误差限为。
7．求方程的两个根，使它至少具有4位有效数字（）。
解：，
故方程的根应为
故
具有5位有效数字
具有5位有效数字
8．当N充分大时，怎样求？
解
设。
则
9．正方形的边长大约为了，应怎样测量才能使其面积误差不超过？
解：正方形的面积函数为
证明：
若
令，可得
当时，
当时，
又，故
得证。
10。证明切xx多项式满足微分方程
证明：
切xx多项式为
从而有
得证。
11。假设在上连续，求的零次最佳一致逼近多项式？
解：
在闭区间上连续
存在，使
取
则和是上的2个轮流为“正”、“负”的偏差点。
由切xx定理知
P为的零次最佳一致逼近多项式。
12。选取常数，使达到极小，又问这个解是否唯一？
解：
令
则在上为奇函数
又的最高次项系数为1，且为3次多项式。
与0的偏差最小。
从而有
13。求在上的最佳一次逼近多项式，并估计误差。
解：
于是得的最佳一次逼近多项式为
清华大学数值分析实验报告
数值分析实验报告一、 实验3。
1 题目：考虑线性方程组b Ax =，n n R A ⨯∈，n R b ∈,编制一个能自动选取主元,又能手动选取主元的求解线性代数方程组的Gauss 消去过程。
（1）取矩阵⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎣⎡=6816816816 A ，⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎣⎡=1415157 b ,则方程有解()T x 1,,1,1*⋯=。
取10=n 计算矩阵的条件数。
分别用顺序Gauss 消元、列主元Gauss 消元和完全选主元Gauss 消元方法求解，结果如何？（2）现选择程序中手动选取主元的功能，每步消去过程都选取模最小或按模尽可能小的元素作为主元进行消元，观察并记录计算结果，若每步消去过程总选取按模最大的元素作为主元,结果又如何？分析实验的结果。
（3）取矩阵阶数n=20或者更大,重复上述实验过程,观察记录并分析不同的问题及消去过程中选择不同的主元时计算结果的差异，说明主元素的选取在消去过程中的作用.(4）选取其他你感兴趣的问题或者随机生成的矩阵,计算其条件数，重复上述实验，观察记录并分析实验的结果。
1. 算法介绍首先,分析各种算法消去过程的计算公式， 顺序高斯消去法:第k 步消去中，设增广矩阵B 中的元素()0k kk a ≠(若等于零则可以判定系数矩阵为奇异矩阵，停止计算），则对k 行以下各行计算()(),1,2,,k ikik k kka l i k k n a ==++,分别用ik l -乘以增广矩阵B 的第k 行并加到第1,2,,k k n ++行，则可将增广矩阵B 中第k 列中()k kka 以下的元素消为零；重复此方法，从第1步进行到第n-1步,则可以得到最终的增广矩阵，即()()(),n n n B Ab ⎡⎤=⎣⎦； 列主元高斯消去法：第k 步消去中，在增广矩阵B 中的子方阵()()()()k kkkknk k nknn a a a a ⎡⎤⎢⎥⎢⎥⎢⎥⎣⎦中,选取()k k i k a 使得()(k)max k k i k ik k i na a ≤≤=，当k i k ≠时，对B 中第k 行与第k i 行交换,然后按照和顺序消去法相同的步骤进行。
数值分析第一次大作业(宋宝瑞)
1、猜想ln(cond(Hn))与n之间在n较小时满足线性关系，在n较大时则不稳定，因此，设n分别为10,20,50,100，如下图所示，可以看出在n<14时，ln(cond(Hn))呈线性稳定增长，n>14时，函数曲线开始趋于平稳并不断波动，波动幅度较小，并有上升趋势。
程序：N=100;y=zeros(N,1);for n=1:NHn=hilb(n);condn=cond(Hn);y(n)=log(condn);endx=1:N;subplot(224);plot(x,y);xlabel('n');ylabel('ln(cond(Hn))');title('ln(cond(Hn))--n(n=100)');grid ;2、可以看出，log(cond(Hnpre)/cond(Hn))在y=0上下波动。
对比整个曲线图，可以发现，经过预处理，条件数有所下降，但下降值不大。
程序：x=1:100;y=zeros(100,1);y1=zeros(100,1);y2=zeros(100,1);for n=1:100Hn=hilb(n);inv_D=zeros(n);for i=1:ninv_D(i,i)=1/sqrt(Hn(i,i));endHnpre=inv_D*Hn*inv_D;y(n)=log(cond(Hnpre)/cond(Hn));y1(n)=log(cond(Hnpre));y2(n)=log(cond(Hn));end%plot(x,y1);plot(x,y,'r',x,y1,'g',x,y2,'k');legend('y','y1','y2');xlabel('n');ylabel('log(cond(Hnpre)/cond(Hn))');title('log(cond(Hnpre)/cond(Hn))--n(n=100)');grid on;3、设b所有值为1，n x1 x2 x34 -4.0000E+00 -4.0000E+00 -4.0000E+00 6.0000E+01 6.0000E+01 6.0000E+01 -1.8000E+02 -1.8000E+02 -1.8000E+02 1.4000E+02 1.4000E+02 1.4000E+0255.0000E+00 5.0000E+00 5.0000E+00 -1.2000E+02 -1.2000E+02 -1.2000E+026.3000E+02 6.3000E+02 6.3000E+02 -1.1200E+03 -1.1200E+03 -1.1200E+03 6.3000E+02 6.3000E+02 6.3000E+026 -6.0000E+00 -6.0000E+00 -6.0000E+00 2.1000E+02 2.1000E+02 2.1000E+02 -1.6800E+03 -1.6800E+03 -1.6800E+03 5.0400E+03 5.0400E+03 5.0400E+03 -6.3000E+03 -6.3000E+03 -6.3000E+03 2.7720E+03 2.7720E+03 2.7720E+0377.0000E+00 7.0000E+00 7.0000E+00 -3.3600E+02 -3.3600E+02 -3.3600E+02 3.7800E+03 3.7800E+03 3.7800E+03-1.6800E+04 -1.6800E+04 -1.6800E+04 3.4650E+04 3.4650E+04 3.4650E+04 -3.3264E+04 -3.3264E+04 -3.3264E+04 1.2012E+04 1.2012E+04 1.2012E+048 -8.0000E+00 -8.0000E+00 -8.0000E+00 5.0400E+02 5.0400E+02 5.0400E+02 -7.5600E+03 -7.5600E+03 -7.5600E+03 4.6200E+04 4.6200E+04 4.6200E+04 -1.3860E+05 -1.3860E+05 -1.3860E+05 2.1622E+05 2.1622E+05 2.1622E+05 -1.6817E+05 -1.6817E+05 -1.6817E+05 5.1480E+04 5.1480E+04 5.1480E+0499.0000E+00 9.0000E+00 9.0001E+00 -7.2000E+02 -7.2000E+02 -7.2000E+02 1.3860E+04 1.3860E+04 1.3860E+04 -1.1088E+05 -1.1088E+05 -1.1088E+05 4.5045E+05 4.5045E+05 4.5045E+05 -1.0090E+06 -1.0090E+06 -1.0090E+06 1.2613E+06 1.2613E+06 1.2613E+06 -8.2368E+05 -8.2368E+05 -8.2368E+05 2.1879E+05 2.1879E+05 2.1879E+0510 -9.9980E+00 -9.9987E+00 -1.0001E+01 9.8983E+02 9.8989E+02 9.9005E+02 -2.3756E+04 -2.3758E+04 -2.3761E+04 2.4021E+05 2.4022E+05 2.4025E+05 -1.2611E+06 -1.2612E+06 -1.2613E+06 3.7833E+06 3.7835E+06 3.7839E+06 -6.7260E+06 -6.7263E+06 -6.7269E+06 7.0006E+06 7.0008E+06 7.0015E+06 -3.9378E+06 -3.9380E+06 -3.9383E+06 9.2370E+05 9.2373E+05 9.2380E+05111.0927E+01 1.0988E+01 1.0992E+01 -1.3124E+03 -1.3188E+03 -1.3192E+03 3.8410E+04 3.8580E+04 3.8587E+04 -4.7823E+05 -4.8015E+05 -4.8022E+05 3.1396E+06 3.1513E+06 3.1515E+06 -1.2060E+07 -1.2102E+07 -1.2102E+072.8484E+07 2.8575E+07 2.8576E+07 -4.1864E+07 -4.1989E+07 -4.1990E+073.7293E+07 3.7398E+07 3.7398E+07 -1.8420E+07 -1.8469E+07 -1.8469E+07 3.8689E+06 3.8786E+06 3.8785E+0612 -9.6349E+00 -1.2768E+01 -1.1739E+01 1.4311E+03 1.8176E+03 1.6833E+03-5.1063E+04 -6.3363E+04 -5.9039E+047.7783E+05 9.4708E+05 8.8711E+05-6.3022E+06 -7.5528E+06 -7.1069E+063.0320E+07 3.5851E+07 3.3868E+07-9.1788E+07 -1.0728E+08 -1.0171E+081.7935E+082.0753E+08 1.9736E+08-2.2572E+08 -2.5889E+08 -2.4688E+081.7662E+082.0099E+08 1.9214E+08-7.8125E+07 -8.8290E+07 -8.4591E+071.4921E+07 1.6757E+07 1.6087E+07可以看出，当n<9时，三种方法计算出的结果几乎一致，当n>9时，三种方法计算出的结果相差越来越大。
清华大学高等数值分析(李津)所有作业答案合集
20130917题目求证：在矩阵的LU 分解中，111n n Tn ij i j j i j L I e e α-==+⎛⎫=- ⎪⎝⎭∑∑证明：在高斯消去过程中，假设0jj a ≠ ，若a=0，可以通过列变换使得前面的条件成立，这里不考虑这种情况。
对矩阵A 进行LU 分解，()()()()()1111111L M n M M M n ---=-=∙∙-………… ，其中()1n Tn ij i j i j M j I e e α=+⎛⎫=+ ⎪⎝⎭∑ ，i e 、j e 为n 维线性空间的自然基。
()M j 是通过对单位阵进行初等变换得到，通过逆向的变换则可以得到单位阵，由此很容易得到()M j 的逆矩阵为1n Tn ij i j i j I e e α=+⎛⎫- ⎪⎝⎭∑。
故111n n T n ij i j n j i j L I e e I α-==+⎛⎫⎛⎫=- ⎪ ⎪ ⎪⎝⎭⎝⎭∏∑上式中的每一项均是初等变换，从右向左乘，则每乘一次相当于对右边的矩阵进行一次向下乘法叠加的初等变换。
由于最初的矩阵为单位阵，变换从右向左展开，因而每一次变换不改变已经更新的数据，既该变换是从右向左一列一列更新数据，故11nn Tn ij i j j i j L I e e α==+⎛⎫=- ⎪⎝⎭∑∑。
数学证明：1nTi j i ji j ee α=+⎛⎫ ⎪⎝⎭∑具有,000n j jA -⎛⎫ ⎪⎝⎭ 和1,1000n j n j B -+-+⎛⎫⎪⎝⎭ 的形式，且有+1,-11,10000=000n j j n j n j AB --+-+⎛⎫⎛⎫⎪⎪⎝⎭⎝⎭ 而11n n T ij i j j k i j e e α-==+⎛⎫ ⎪⎝⎭∑∑具有1,1000n k n k B -+-+⎛⎫⎪⎝⎭的形式，因此：1311111211121==n n n n n n T T T n ij i j n ij i j n ik i k j i j j i j k n i k n n T n i i n ik i i i k L I e e I e e I e e I e e I e ααααα---==+==+=-=+==+⎡⎤⎛⎫⎛⎫⎛⎫⎛⎫⎛⎫⎛⎫=---⎢⎥ ⎪ ⎪ ⎪ ⎪ ⎪⎪ ⎪ ⎪⎝⎭⎝⎭⎢⎥⎝⎭⎝⎭⎝⎭⎝⎭⎣⎦⎛⎫⎛⎫⎛⎫=-- ⎪ ⎪ ⎝⎭⎝⎝⎭∏∑∏∑∑∑∑∑……11211n n n T Tk n ik i kk k i k e I e e α--===+⎛⎫⎛⎫=- ⎪⎪ ⎪⎭⎝⎭⎝⎭∑∑∑#20130924题目一问：能否用逐次householder 相似变换变实矩阵A 为上三角矩阵，为什么？解：不能用逐次householder 相似变换变A 为上三角矩阵，原因如下：A 记作：()12=,,n A a a a ……， ，存在householder 阵1H s.t. 1111H a e α= ，则()()()111111111111111111111,,,0T Th H AH H a A H e H A H e H A H h H A H ααα⎛⎫'''=== ⎪⎪'⎝⎭⎛⎫''=+ ⎪ ⎪⎝⎭11H A H ''第一列的元素不能保证为1e 的倍数，故无法通过householder 变换实现上三角化。
数值分析第五版习题答案清华大学出版社
1．设 x 0 , x 的相对误差为 ，求 ln x 的误差。
解：近似值
*
x 的相对误差为
* e* x * x
= er
x*
x*
而 ln x 的误差为 e ln x *
ln x * ln x
1 e*
x*
进而有 (ln x*)
2．设 x 的相对误差为 2%，求 xn 的相对误差。
解：设 f ( x )
cos x 的近似
值时，采用的线性插值法插值余项不为 0，也会有一定的误差。因此，总误差界的计算应综
合以上两方面的因素。
当 0 x 90 时，
令 f ( x ) cos x
1 取 x0 0, h ( )
60
1 60 180
令 xi x 0 ih , i 0,1,..., 5400
10800
则 x5400
3
2．给出 f ( x ) ln x 的数值表
X
0.4
lnx
-0.916291
用线性插值及二次插值计算
解：由表格知，
0.5
0.6
-0.693147
-0.510826
ln 0.54 的近似值。
x0 0.4, x1 0.5, x 2 0.6, x3 0.7, x4 0.8; f ( x 0 ) 0.916291, f ( x1 ) 0.693147 f ( x 2 ) 0.510826, f ( x3 ) 0.356675 f ( x 4 ) 0.223144
2 1.41
( y 0 *)
1
2
10
2
又 y n 10 yn 1 1
y1 10 y0 1
( y1*) 10 ( y 0 *) 又 y 2 10 y1 1
清华大学数值分析A往年试题回顾
5 道大题，若干小题，卷面成绩满分70
1.(1)求f(x)=sqrt(1-x A2)在span{1,x,xA2}上，权函数为rou=1/sqrt(1-x A2)的最佳平方逼近多项式
⑵求证高斯型求积公式中的A(k)满足A(k)= / p(x)l(x)dx= / p(x)lA其(X)dXk)为Lagrange多项
式
2.(1)Ax=b中A非奇异，则用J法、GS法、SOR法、SSOR法求解等价方程ATAx=ATb各种方法的收敛性怎样？(其中0<w<2)
(2)A严格对角占优，求证其有唯一的LU分解，对称矩阵[3 1 0;1 3 1;0 1 3]求其cholysky分解
3.(1)写出用Lanczos方法计算某矩阵第一列的a和B
⑵已知矩阵[3 0 0;0 3 2;0 2 3]，求其QR分解，计算一步H'=RQ
4(1)f(x)=[x2A2-x1A2-x1 其精确解为x*=[0 0 0],写出牛顿法的计算公式sin(x1A2)-x2];
(2)已知G(x)=[x2A2-x1A2 sin(x1A2)];
给出区域D 使得在此区域内的初始值可以收敛到精确解，并说明原因
5.(1)线性2 步法-0.5y(n)-0.5y(n+1)+y(n+2)=h/2*(f(n)+f(n+1)+f(n+2))，计算其局部阶段误差的阶数若h=0.1,判断其稳定性
⑵已知R(z)的稳定函数是exp(z)的pade(1,2)逼近多项式，计算其稳定域，是否是A-稳定？(pade 逼近的计算公式卷子上给了)。
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
7、设y0=28，按递推公式
y n=y n−1−
1
100
√783，n=1,2,…
计算y100，若取√783≈27.982，试问计算y100将有多大误差？
答：y100=y99−1
100√783=y98−2
100
√783=⋯=y0−100
100
√783=28−√783
若取√783≈27.982，则y100≈28−27.982=0.018，只有2位有效数字，y100的最大误差位0.001
10、设f(x)=ln(x−√x2−1)，它等价于f(x)=−ln(x+√x2−1)。
分别计算f(30)，开方和对数取6位有效数字。
试问哪一个公式计算结果可靠？为什么？
答：√x2−1≈29.9833
则对于f(x)=ln(x−√x2−1)，f(30)≈−4.09235
对于f(x)=−ln(x+√x2−1)，f(30)≈−4.09407
而f(30)=ln(30−√302−1)，约为−4.09407，则f(x)=−ln(x+√x2−1)计算结果更可靠。
这是因为在公式f(x)=ln(x−√x2−1)中，存在两相近数相减（x−√x2−1）的情况，导致算法数值不稳定。
11、求方程x2+62x+1=0的两个根，使它们具有四位有效数字。
答：x12=−62±√622−4
2
=−31±√312−1
则
x1=−31−√312−1≈−31−30.98=−61.98
x2=−31+√312−1=
1
31+√312−1
≈−
1
31+30.98
≈−0.01613
12.（1）、计算√101.1−√101，要求具有4位有效数字
答：√101.1−√101=
√101.1+√101≈0.1
10.05+10.05
≈0.004975
14、试导出计算积分I n=∫x n
4x+1dx
1
的一个递推公式，并讨论所得公式是否计算稳定。
答：I n=∫x n
4x+1dx
1 0=∫
1
4
(4x+1)x n−1−1
4
x n−1
4x+1
dx=
1
1
4
∫x n−1
1
dx−1
4
∫x n−1
4x+1
dx
1
=
1
4n
−
1
4
I n−1，n=1，2…
I0=∫
1
4x+1
dx=
ln5
4
1
记εn为I n的误差，则由递推公式可得
εn=−1
4
εn−1=⋯=(−
1
4
)nε0
当n增大时，εn是减小的，故递推公式是计算稳定的。