下载文档原格式
数 值 分 析（B ） 大 作 业（一）姓名： 学号： 电话：1、算法设计：①求1λ、501λ和s λ的值：s λ：s λ表示矩阵的按模最小特征值，为求得s λ直接对待求矩阵A 应用反幂法即可。
1λ、501λ：若矩阵A 的特征值满足关系 1n λλ<<且1n λλ≠，要求1λ、及501λ时，可按如下方法求解： a ． 对矩阵A 用幂法，求得按模最大的特征值1m λ。
b ． 按平移量1m λ对矩阵A 进行原点平移得矩阵1m BA I λ=+，对矩阵B 用反幂法求得B 的按模最小特征值2m λ。
c ． 321m m m λλλ=-则：113min(,)m m λλλ=，13max(,)n m m λλλ=即为所求。
②求和A 的与数5011140k k λλμλ-=+最接近的特征值ik λ（k=0，1，…39）：求矩阵A 的特征值中与P 最接近的特征值的大小，采用原点平移的方法：先求矩阵 B=A-PI 对应的按模最小特征值k β，则k β+P 即为矩阵A 与P 最接近的特征值。
在本次计算实习中则是先求平移矩阵k B A I μ=-，对该矩阵应用反幂法求得s λ，则与k μ最接近的A 的特征值为：s P λ+重复以上过程39次即可求得ik λ（k=0，1，…39）的值。
③求A 的（谱范数）条件数2cond()A 和行列式det A ：在（1）中用反幂法求矩阵A 的按模最小特征值时，要用到Doolittle 分解方法，在Doolittle 分解完成后得到的两个矩阵分别为L 和U ，则A 的行列式可由U 阵求出，即：det(A)=det(U)。
求得det(A)不为0，因此A 为非奇异的实对称矩阵，则： max 2()scond A λλ=，max λ和s λ分别为模最大特征值与模最小特征值。
2、程序源代码：#include "Stdio.h"#include "Conio.h"#include "math.h"//****************************************************************************// // 在存储带状矩阵时，下面的几个量在程序中反复用到，为方便编程故把它们定义成宏.// // M ：转换后的矩阵的行数，M=R+S+1。
数值分析第一次作业及参考答案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、幂法求解λ1、λ501幂法主要用于计算矩阵的按模最大的特征值和相应的特征向量，即对于|λ1|≥|λ2|≥.....≥|λn|可以采用幂法直接求出λ1，但在本题中λ1≤λ2≤……≤λ501，我们无法判断按模最大的特征值。
但是由矩阵A的特征值条件可知|λ1|和|λ501|之间必然有一个是最大的，通过对矩阵A使用幂法迭代一定次数后得到满足精度ε=10−12的特征值λ0，然后在对矩阵A做如下的平移：B=A-λ0I由线性代数(A-PI)x=(λ-p)x可得矩阵B的特征值为：λ1-λ0、λ2-λ0…….λ501-λ0。
对B矩阵采用幂法求出B矩阵按模最大的特征值为λ∗=λ501-λ0，所以λ501=λ∗+λ0，比较λ0与λ501的大小，若λ0>λ501则λ1=λ501，λ501=λ0；若λ0<λ501，则令t=λ501，λ1=λ0，λ501=t。
求矩阵M按模最大的特征值λ的具体算法如下：任取非零向量u0∈R nηk−1=u T(k−1)∗u k−1y k−1=u k−1ηk−1u k=Ay k−1βk=y Tk−1u k(k=1,2,3……)当|βk−βk−1||βk|≤ε=10−12时，迭终终止，并且令λ1=βk2、反幂法计算λs和λik由已知条件可知λs是矩阵A 按模最小的特征值，可以应用反幂法直接求解出λs。
使用带偏移量的反幂法求解λik，其中偏移量为μk=λ1+kλ501−λ140(k=1,2,3…39)，构造矩阵C=A-μk I，矩阵C的特征值为λik−μk，对矩阵C使用反幂法求得按模最小特征值λ0，则有λik=1λ0+μk。
求解矩阵M按模最小特征值的具体算法如下：任取非零向量u 0∈R n ηk−1= u T (k−1)∗u k−1y k−1=u k−1ηk−1 Au k =y k−1βk =y T k−1u k (k=1,2,3……)在反幂法中每一次迭代都要求解线性方程组Au k =y k−1，当K 足够大时，取λn =1βk 。
问题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代表端点的一阶导。
MATLAB作业1. 判断如下命题是否正确(a) 一个问题的病态性如何，与求解它的算法有关系。
(错) (b) 无论问题是否病态，好的算法都会得到它好的近似解。
（错） (c) 计算中使用更高的精度，可以改善问题的病态性。
(错) (d) 用一个稳定的算法计算一个良态的问题，一定会得到它好的近似解。
(错) (e) 浮点数在整个数轴上是均匀分布的 (错) (f) 浮点数的加法满足结合律 (错) (g) 浮点数的加法满足交换律 (错) (h) 浮点数构成有限集合 (对） (i) 用一个收敛的算法计算一个良态的问题，一定会得到它好的近似解 (错)2. 函数sinx 有幂级数展开利用幂级数计算sinx 的Matlab 程序为 function s = powersin(x)% POWERSIN. Power series for sin(x).% POWERSIN(x) tries to compute sin(x) from a power series s = 0; t = x; n = 1;while s + t ~= s; s = s + t;t = -x^2/((n+1)*(n+2))*t; n = n + 2; end(a) 解释上述程序的终止准则；(b) 对于，计算的精度是多少？分别需要计算多少项？答：（a ）当t 小于计算机的计算精度时，上述程序将终止。
（b ） x=/2π; n=23； s=1.0000x=11/2π; n=75； s= -1.0000x =21/2π; n=121； s= 0.9999/2,11/2,21/2x πππ=3. 考虑数列 ，它的统计平均值定义为它的标准差数学上等价于作为标准差的两种算法，你如何评价它们的得与失？第一种算法共进行了n 次乘方运算，2n 次求和运算，第二种算法进行了2n 次乘方和2n 次求和，运算次数较多。
而且第二种算法中2i x 与误差较为接近，易造成舍入2x n 。
习题（一）1. 指出四舍五入得到的下列各数有几位有效数字：x 1∗=7.8673,x 2∗=8.0916,x 3∗=0.06213,x 4∗=0.07800,x 5∗=90×103,x 6∗=2.0×10−4解：由有效数字定义得：x 1∗，x 2∗具有5位有效数字x 3∗，x 4∗具有4位有效数字x 5∗，x 6∗具有2位有效数字.2. 设准确值为x=3.78695，y=10,它们的近似值分别为x 1∗=3.7869，x 2∗=3.7870及y 1∗=9.9999，y 2∗=10.1，y 3∗=10.0001，试分析x 1∗，x 2∗，y 1∗，y 2∗，y 3∗分别具有几位有效数字. 解：x 1∗=3.7869=x 1∗=0.37869×101，k 1=1|x 1∗−x|=|3.7869−3.78695|=0.00005≤0.5×10−4=0.5×101−5, 即x 1∗具有5位有效数字；同理，x 2∗=3.7870=0.37870×101，k 2=1|x 2∗−x|=|3.7870−3.78695|=0.00005≤0.5×101−5，所以x 2∗具有5位有效数字； 将y 1∗，y 2∗，y 3∗分别写成y=±10k ×0.α1α2...αn 的表示形式，有：y 1∗=9.9999=0.99999×101，k 3=1；y 2∗=10.1=y 2∗=0.101×102，k 4=2；y 3∗=10.0001=0.100001×102，k 5=2；|y 1∗−y |=|9.9999−10|=0.0001=0.1×10−3≤0.5×101−4,n=4;|y 2∗−y |=|10.1−10|=0.1≤0.5×102−2,n=2;|y 3∗−y |=|10.0001−10|=0.0001=0.1×10−3≤0.5×102−5,n=5;所以y 1∗，y 2∗，y 3∗分别具有4，2，5位有效数字.8.为了使√11的近似值的相对误差不超过0.1%，问至少应取几位有效数字. 解：√11=0.3316624…=0.α1α2...αn ×10k ，α1=3，设x ∗有n 位有效数字，又因为|E x ∗|比值比较小， 故可用E r ∗(x ∗)= |E(x ∗)x ∗|代替相对误差E r ∗(x ∗)，用εr ∗=εx ∗代替相对误差限εr 所以εr ∗≤12α1×10−n+1=16×10−n+1 令16×10−n+1≤0.1%，解得n ≥3.22即至少应取4位有效数字.12.如何计算下列函数值才比较精确.(1)11+2x −11+x ，对|x|≪1; (2)√x +1x −√x −1x ，对x ≫1;(3)∫dx 1+x 2N+1N，其中N 充分大； (4)1−cos xsin x ，对|x|≪1;(5)ln(30−√302−1)(开平方用6位函数表)；解:(1)原式=1+x−(1+2x)(1+2x)(1+x)=−x (1+2x)(1+x)； (2)原式=x+1x −(x−1x )√x+1x +√x−1x =2x √x+1x +√x−1x ；(3)原式=arc tan x|NN+1=arc tan N +1−arc tan N =arc tan N+1−N 1+N(N+1)=arc tan 11+N(N+1)； (4)原式=2sinx 222sin x 2cos x 2=tan x2; (5)原式=30+√302−1=−ln(30+√302−1)令f(x)=ln(x −√x 2−1)，则f(30)=ln(30−√302−1)=ln(30−√899),记a=30−√899 若用6位开方函数表，则有a ∗=30−29.9833=0.0167,故有ε(a ∗)=0.5×10−4, 而f(30)≈ln a ∗,于是ε(f (30))=ε(ln a ∗)≈|1a ∗|ε(a ∗)=0.50.0167×10−4≈0.003; 又因为f(x)等价于f(x)=-ln(x +√x 2−1)，则f (30)=-ln(30+√899)，记b=30+√899 同理b ∗=59.9833,进而ε(b ∗)=(2×10−4)−1,对f (30)≈ln b ∗ε(f (30))=ε(ln b ∗)≈|1b ∗|ε(b ∗)=0.559.9833×10−4≈0.834×10−6。
数值分析期末考试一、 设80~=x ，若要确保其近似数的相对误差限为0.1%，则它的近似数x 至少取几位有效数字？（4分）解：设x 有n 位有效数字。
因为98180648=<<=，所以可得x 的第一位有效数字为8（1分） 又因为21101011000110821--⨯=<⨯⨯≤n ε，令321=⇒-=-n n ，可知x 至少具有3位有效数字（3分）。
二、求矩阵A 的条件数1)(A Cond （4分）。
其中⎥⎦⎤⎢⎣⎡=4231A 解：⎥⎦⎤⎢⎣⎡--=-5.05.1121A （1分） 1A =7（1分） 2711=-A （1分）249)(1=A Cond （1分）三、用列主元Gauss 消元法法求解以下方程组（6分）942822032321321321=++-=++--=+-x x x x x x x x x解：→⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡---→⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡----→⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡----5.245.2405.35.230914220321821191429142821120321 ⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡---→⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡---8175835005,245.24091425.33.2305.245.2409142(4分) 等价三角方程组为：⎪⎪⎩⎪⎪⎨⎧-=-=+-=++,8175835,5.245.24,942332321x x x x x x （1分）回代得1,3,5123==-=x x x （1分）四、设.0,2,3,1,103)(3210234=-===-+-=x x x x x x x x f 1）求以3210,,,x x x x 为节3次Lagrange 多项式；（6分） 2）求以3210,,,x x x x 为节3次Newton 多项式；（6分）3）给出以上插值多项式的插值余项的表达式（3分）解：由0,2,3,13210=-===x x x x 可得10)(,34)(,1)(,11)(3210-==-=-=x f x f x f x f即得： +------+------=))()(())()(()())()(())()(()()(312101320130201032103x x x x x x x x x x x x x f x x x x x x x x x x x x x f x L=------+------))()(())()(()())()(())()(()(23130321033212023102x x x x x x x x x x x x x f x x x x x x x x x x x x x f+-+--+-⨯-+-+--+-⨯-)03)(23)(13()0)(2)(1()1()01)(21)(31()0)(2)(3(11x x x x x x326610.)20)(30)(10()2)(3)(1()10()02)(32)(12()0)(3)(1(34x x x x x x x x x -+--=+--+--⨯-+---------⨯2）计算差商表如下：i x )(i x f 一阶差商 二阶差商 三阶差商1 -11 3 -1 5 -2 34 -7 4 0-10-225-1则=+-----+-+-=)2)(3)(1()3)(1(4)1(511)(3x x x x x x x N326610x x x -+--3）)2)(3)(1())()()((!4)()(3210)4(3+--=----=x x x x x x x x x x x x f x R ξ五、给定方程组b Ax =，其中⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡=100131w w w w A 。
北京航空航天大学数值分析大作业一学院名称自动化专业方向控制工程学号ZY*******学生姓名许阳教师孙玉泉日期2021 年11月26 日设有501501⨯的实对称矩阵A ，⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎣⎡=5011A a b c b c c b c b a其中，064.0,16.0),501,,2,1(64.0)2.0sin()024.064.1(1.0-==⋅⋅⋅=--=c b i e i i a ii 。
矩阵A 的特征值为)501,,2,1(⋅⋅⋅=i i λ，并且有||min ||,501150121i i s λλλλλ≤≤=≤⋅⋅⋅≤≤1λ，501λ和s λ的值。
A 的与数4015011λλλμ-+=kk 最接近的特征值)39,,2,1(⋅⋅⋅=k k i λ。
A 的(谱范数)条件数2)A (cond 和行列式detA 。
一 方案设计1 求1λ，501λ和s λ的值。
s λ为按模最小特征值，||min ||5011i i s λλ≤≤=。
可使用反幂法求得。
1λ，501λ分别为最大特征值及最小特征值。
可使用幂法求出按模最大特征值，如结果为正，即为501λ，结果为负，那么为1λ。
使用位移的方式求得另一特征值即可。
2 求A 的与数4015011λλλμ-+=kk 最接近的特征值)39,...,2,1(=k k i λ。
题目可看成求以k μ为偏移量后，按模最小的特征值。
即以k μ为偏移量做位移，使用反幂法求出按模最小特征值后，加上k μ，即为所求。
3 求A 的(谱范数)条件数2)(A cond 和行列式detA 。
矩阵A 为非奇异对称矩阵，可知，||)(min max2λλ=A cond(1-1)其中m ax λ为按模最大特征值，min λ为按模最小特征值。
detA 可由LU 分解得到。
因LU 均为三角阵，那么其主对角线乘积即为A 的行列式。
二 算法实现1 幂法使用如下迭代格式：⎪⎪⎩⎪⎪⎨⎧⋅===⋅⋅⋅=------||max |)|sgn(max ||max /),,(111111)0()0(10k k k k k k k k Tn u u Ay u u u y u u u β任取非零向量 (2-1)终止迭代的控制理论使用εβββ≤--||/||1k k k ， 实际使用εβββ≤--||/||||||1k k k(2-2)由于不保存A 矩阵中的零元素，只保存主对角元素a[501]及b,c 值。
第一章 误差与算法1. 误差分为有__模型误差___, _观测误差___, __方法误差____, ___舍入误差____, Taylor 展开式近似表达函数产生的误差是_方法误差 .2. 插值余项是插值多项式的 方法误差。
0.2499作为1/4的近似值, 有几位有效数字？00.24990.249910,0m =⨯=即，031|0.2499|0.00010.5100.510,34m n n ---=<⨯=⨯=即22 3.1428751...,7=作为圆周率的近似值，误差和误差限分别是多少，有几位有效数字？2133.142875 3.14159260.00126450.5100.510---=<⨯=⨯有3位有效数字.* 有效数字与相对误差的关系3. 利用递推公式计算积分110,1,2,...,9n x n I x e dx n -==⎰错误!未找到引用源。
, 建立稳定的数值算法。
该算法是不稳定的。
因为:11()()...(1)!()n n n I n I n I εεε-=-==-111n n I I n n -=-， 10110I =4. 衡量算法优劣的指标有__时间复杂度,__空间复杂度_.时间复杂度是指: , 两个n 阶矩阵相乘的乘法次数是 , 则称两个n 阶矩阵相乘这一问题的时间复杂度为 .二 代数插值1.根据下表数据建立不超过二次的Lagrange 和Newton 插值多项式, 并写出误差估计式, 以及验证插值多项式的唯一性。
x 0 1 4f(x) 1 9 3Lagrange:设0120120,1,4;()1()9()3x x x f x f x f x ======则，， 对应 的标准基函数 为:1200102()()(1)(x 4)1()(1)(x 4)()()(01)(04)4x x x x x l x x x x x x ----===------ 1()...l x =2()...l x =因此, 所求插值多项式为:220()()()....i i i P x f x l x ===∑ (3)2()()(0)(1)(x 4)3!f R x x x ξ=--- Newton:构造出插商表:xi f(xi ) 一 二 三0 11 9 84 3 -2 -5/2所以, 所求插值多项式为:2001001201()()[,]()[,,]()()518(0)(0)(1)2...P x f x f x x x x f x x x x x x x x x x =+-+--=+----=插值余项: 2()[0,1,4,](0)(1)(x 4)R x f x x x =---2. 已知函数f(0)=1,f(1)=3,f(2)=7,则f[0,1]=___2________, f[0,1,2]=____1______)('],[000x f x x f =3.过0,1两节点构造三次Hermite 插值多项式, 使得满足插值条件: f(0)=1. .’(0)=... f(1.=2. .’(1)=1设0101010,1,()1()2'()0,'()1x x f x f x f x f x ======则，， 写出插商表:xi f(xi) 一 二 三0 10 1 01 a 1 11 a 1 0 a-1因此, 所求插值多项式为:插值余项:222()[0,0,1,1,](1)R x f x x x =-4.求f(x)=sinx 在[a,b]区间上的分段线性插值多项式, 并写出误差估计式。
数值分析课后作业：习题一1.在字长为3的十进制计算机上计算f （3.33）和g （3.33），其中f(x)=x 4-x 3+3x 2+x-2,g(x)=(((x-1)x+3)x+1)x-2解： m=3; f=@(x)digit(digit(x^4,m)- digit(x^3,m)+ digit(3*x^2,m)+ digit(x-2,m),m); g=@(x)digit(digit(digit( digit(digit(digit( (x-1)*x,m)+3,m)*x,m)+1,m)*x,m)-2,m); f(3.33) g(3.33) 有ans = 121 ans =121 2.下列各近似值的绝对误差限都是1021⨯-3，试指出它们各有几位有效数字：x=1.00052， y=0.05， z=0.00052.解：当 x=1.00052时， 由丨X*—X 丨 ≤0.5×10-3 得 x=1.00052 有四位有效数字； 同理 y=0052 有两位有效数字 Z=0.00052有零位有效数字 3,计算圆的面积，要使其相对误差限为1%，问测量半径r 允许的相对误差限是多少？ 解：设圆的面积为S ， 由题意有|e(S)|≤1%。
又S=πr 2 dS=2πr dr 所以 dS/S=(2πrdr)/(πr 2)=2(dr/r)∴|e(r)|≈21|e(S)|≤0.5×1%=0.5% 11.数组与矩阵是Matlab 编程的基础，试学习Matlab 的数组与矩阵的表示方法，并举例介绍数组、矩阵的常见运算. 解：>> syms a b c d; >> a=[1 2 3];>> b=[4 5 6];>> a+bans =5 7 9>> b-aans =3 3 3>> a.*bans =4 10 18 >> a.^2 ans = 1 4 9>> c=[1 2 3;1 2 3;1 2 3];>> d=[4 5 6;4 5 6;4 5 6];>> cc = 1 2 3 1 2 3 1 2 3d = 4 5 6 4 5 6 4 5 6 >> c+dans =5 7 9 5 7 9 5 7 9>> d-cans = 3 3 33 3 33 3 3 12.学习使用Matlab 命令help 和doc 学习自己感兴趣的Matlab 的运算、函数或命令的用法，并对于任意给定的实数a,b,c,编写Matlab 程序求方程ax 2+bx+c=0的根. 解：x 1=a ac b b b 24)sgn(2---, x 2=1ax c1 x>0 其中 sgn = 0 x=0 -1 x<0 disp('Please input the coefficients of');disp('quadratic equation ax^2+bx+c=0, respectively') a=input('a='); b=input('b='); c=input('c=');m=3; if abs(a)<eps & abs(b)<eps error End if abs(a)<eps disp('Since a=0, quadrtic equation degen erates into a linear equation.') disp('The only solution of the linear equtio n is')x=digit(-c/b,m) return Enddelta=b^2-4*a*c; temp=sqrt(delta); x 1=(-b+temp)/(2*a) ; x 2=(-b-temp)/(2*a) ;err1=abs(a*x 1^2+b*x 1+c) ; err2=abs(a*x 2^2+b*x 2+c) ; if b>0x 1=(-b-temp)/(2*a) End if b<0x 1=(-b+temp)/(2*a) End if b=0x 1=temp/(2*a) Endx 2=c/(a*x 1)err1=abs(a*x 1^2+b*x 1+c) err2=abs(a*x 2^2+b*x 2+c) if abs(a)<epsdisp('Since a=0, quadrtic equation degen erates into a linear equation.')disp('The only solution of the linear equtio n is')x=digit(-c/b,m) return Enddelta=digit(digit(b^2,m)-digit(4*digit(a*c,m),m),m);temp=digit(sqrt(delta),m);x 1=digit(digit(-b+temp,m)/digit(2*a,m),m); x 2=digit(digit(-b-temp,m)/digit(2*a,m),m); err1=abs(a*x 1^2+b*x 1+c); err2=abs(a*x 2^2+b*x 2+c); if b>0x 1=digit(digit(-b-temp,m)/digit(2*a,m),m) ; End if b<0x 1=digit(digit(-b+temp,m)/digit(2*a,m),m); End if b=0x 1=digit(temp/digit(2*a,m),m); Endx 2=digit(digit(c/a,m)/x1,m) ; err1=abs(a*x 1^2+b*x 1+c) ; err2=abs(a*x 2^2+b*x 2+c) ; 14分别利用ln (1+x)=11,)1(11≤<--+∞=∑x nx nn n 和ln11...),12...53(2111253<<-++++++=-++x n x x x x x x n ，给出计算ln2的近似方法，编写相应的Matlab 程序，并比较算法运行情况. 解：方法一： x=1; s=0;for k=1:100s=s+(-1)^(k+1)*(x^k)/k; end sq=log(2)err=abs(t-q) ans= t =0.6882 q =0.6931 err = 0.0050方法二x=1/3; s=0;for k=1:2:100 s=s+(x^k)/k; end t=2*s q=log(2)err=abs(t-q) Ans= t =0.6931 q =0.6931 err =2.2204e-16所以方法二较方法一好。
问题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代表端点的一阶导。
《数值分析》计算实习作业《一》北航第一题 设有501501⨯的矩阵⎥⎥⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎢⎢⎣⎡=501500499321a bc b a b cc b a b ccb a bc c b a b c b a A其中.064.0,16.0);501,2,1(64.0)2.0sin()024.064.1(1.0-===--=c b i e i i a i i 矩阵的特征值)501,,2,1( =i i λ满足||min ||,501150121i i s λλλλλ≤≤=<<<试求1. 5011，λλ和s λ的值2. 的与数4015011λλκλμ-+=k 最接近的特征值)39,,2,1( =K κλi3. 的（谱范数）条件数2)A (cond 和行列式A det 要求1. 算法的设计方案（A 的所有零元素都不能存储）2. 全部源程序（详细注释）。
变量为double ，精度-1210=ε，输出为e 型12位有效数字3. 特征值s 5011，，λλλ和)39,,2,1( =K κλi 以及A cond det ,)A (2的值4. 讨论迭代初始向量的选取对计算结果的影响，并说明原因解答：1. 算法设计对于s λ满足||min ||5011i i s λλ≤≤=，所以s λ是按模最小的特征值，直接运用反幂法可求得。
对于5011，λλ，一个是最大的特征值，一个是最小的特征值，不能确定两者的绝对值是否相等，因此必须首先假设||||5011λλ≠，然后运用幂法，看能否求得一个特征值，如果可以求得一个，证明A 是收敛的，求得的结果是正确的，然后对A 进行带原点平移的幂法，偏移量是前面求得的特征值，可以求得另一个特征值，最后比较这两个特征值，较大的特征值是501λ，较小的特征值就是1λ。
如果在假设的前提下，无法运用幂法求得按模最大的特征值，即此时A 不收敛，则需要将A 进行带原点平移的幂法，平移量可以选取1，再重复上述步骤即可求得两个特征值。
作业1.用如下数值表构造不超过3次的插值多项式2. P55 11题.给出概率积分⎰-=xxdxey 022π的数据表用2次插值计算，试问：(1) 当x = 0.472时，积分值等于多少? (2) 当x 为何值时，积分值等于0.5？ 解：(1) 取x 0 = 0.47， x 1 = 0.48, x 2 = 0.4980.4955530040.04093346-80.1809899240.355496540.51166830.50274980.4937452=+=----⨯+----⨯+----⨯==----+----+----≈)48.049.0)(47.049.0()48.0472.0)(47.0472.0()49.048.0)(47.048.0()49.0472.0)(47.0472.0()49.0472.0)(48.047.0()49.0472.0)(48.0472.0()472.0())(())(())(())(())(())(()472.0(2120210221120121210y Lxx xx xxy x x x x x x y x x x x x x x x x x x x y(2)90.4769359350.05272367-80.4362204360.093439170.50274980.51166830.49374520.51166830.50274980.49374520.49 0.51166830.50274980.49374520.50274980.51166830.49374520.48 0.51166830.49374520.50274980.49374520.51166830.50274980.47=+=----⨯+----⨯+----⨯==----+----+----≈))(()5.0)(5.0())(()5.0)(5.0())(()5.0)(5.0()5.0())(())(())(())(())(())(()5.0(212210221120121210Lyyy y yy xy y yy y y xyyyy yy xy y y y y y x3. 证明方程e x +10x －2＝0在区间[0，1]内有一个根，如果使用二分法求该区间内的根，且误差不超过10－6，试问需要二分区间[0，1]多少次？4. 设x t ＝451.01为准确值，x a ＝451.023为x t 的近似值，试求出x a 有效数字的位数及相对误差 作业答案1.解：N 2(x ) = f (0)+f [0，1](x -0)+ f [0，1，2](x -0) (x -1) 1＋1×(x -0) ＋3×(x -0) (x -1)＝3x 2-2x +1 为求得P 3(x )，根据插值条件知，P 3(x )应具有下面的形式 P 3(x )＝N 2(x )＋k (x -0) (x -1) (x -2)，这样的P 3(x )自然满足：P 3(x i )= f (x i )由P 3’(1 )=3P 3’(1 )= N 2’(1 )+k (1-0) (1-2) =N 2’(1 )－k = 4－k=3∴ k =1∴ P 3(x )＝N 2(x )＋ (x -0) (x -1) (x -2)=x 3+1 3. 证明 令f (x )＝e x +10x －2，∵ f (0)=-1<0，f (1)=e+8> 0∴ f (x )= e x +10x －2 =0在[0，1]有根。
数值分析大作业一、算法设计方案1、矩阵初始化矩阵[]501501⨯=ij a A 的下半带宽r=2,上半带宽s=2，设置矩阵[][]5011++s r C ，在矩阵C 中检索矩阵A 中的带内元素ij a 的方法是：j s j i ij c a ,1++-=。
这样所需要的存储单元数大大减少，从而极大提高了运算效率。
2、利用幂法求出5011λλ，幂法迭代格式：0111111nk k k k kk T k k k u R y u u Ay y u ηηβ------⎧∈⎪⎪=⎪=⎨⎪=⎪⎪=⎩非零向量 当1210/-≤-k k βββ时，迭代终止。
首先对于矩阵A 利用幂法迭代求出一个λ，然后求出矩阵B ，其中I A B λ-=(I 为单位矩阵），对矩阵B 进行幂法迭代，求出λ'，之后令λλλ+'=''，比较的大小与λλ''，大者为501λ，小者为1λ。
3、利用反幂法求出ik s λλ,反幂法迭代格式：0111111nk k k k kk T k k k u R y u Au y y u ηηβ------⎧∈⎪⎪=⎪=⎨⎪=⎪⎪=⎩非零向量 当1210/-≤-k k βββ时，迭代终止，1s k λβ=。
每迭代一次都要求解一次线性方程组1-=k k y Au ，求解过程为：（1）作分解LU A =对于n k ,...,2,1=执行[][]s k n r k k k i c c c c c n s k k k j c cc c k s ks k t k s k r i t t s t i k s k i k s k i js j t k s j r k t t s t k j s j k j s j k <+++=-=++=-=+++----=++-++-++-++----=++-++-++-∑∑);,min(,...,2,1/)(:),min(,...,1,:,1,11),,1max(,1,1,1,11),,1max(,1,1,1（2）求解y Ux b Ly ==,（数组b 先是存放原方程组右端向量，后来存放中间向量y))1,...,2,1(/)(:/:),...,3,2(:,1),min(1.1.11),1max(,1--=-===-=+++-++-+--=++-∑∑n n i c x c b x c b x n i b c b b i s t n s i i t t s t i i i ns n n ti r i t t s t i i i使用反幂法，直接可以求得矩阵按模最小的特征值s λ。
数值分析作业第一题一、 算法设计方案利用带状Dollittle 分解，将A[501][501]转存到数组C[5][501]，以节省存储空间1、计算λ1和λ501首先使用幂法求出矩阵的按模最大的特征值λ0：如果λ0>0，则其必为按模最大值，因此λ501=λ0，然后采用原点平移法，平移量为λ501，使用幂法迭代求出矩阵A -λ501I 的按模最大的特征值，其特征值按从小到大排列应为λ1-λ501、λ2-λ501、……、0。
因此A-λ501I 的按模最大的特征值应为λ1-λ501。
又因为λ501的值已求得，由此可直接求出λ1。
2、计算λSλS 为矩阵A 按模最小的特征值，可以通过反幂法直接求出。
3、计算λikλik 是对矩阵A 进行λik 平移后，再用反幂法求出按模最小的特征值λmin ，λik =λik +λmin 。
4、计算矩阵A 的条件数计算cond （A ）2和行列式det(A)矩阵A 的条件数为n12cond λλ）（ A ，其中λ1和λn 分别是矩阵A 的模最大和最小特征值，直接利用上面求得的结果直接计算。
矩阵A 的行列式可先对矩阵A 进行LU 分解后，det(A)等于U 所有对角线上元素的乘积。
二、源程序：#include<math.h>#include<stdio.h>#include<stdlib.h>#include<iostream.h>#define s 2#define r 2int Max(int v1,int v2);int Min(int v1,int v2);int maxt(int v1,int v2,int v3);void storage(double C[5][501],double b,double c);double mifa(double C[5][501]);void LU(double C[5][501]);double fmifa(double C[5][501]);int Max(int v1,int v2) //求两个数的最大值{ return((v1>v2)?v1:v2);}int Min(int v1,int v2) //求两个数最小值{ return ((v1<v2)?v1:v2);}int maxt(int v1,int v2,int v3) //求三个数最大值{ int t;if(v1>v2) t=v1;else t=v2;if(t<v3) t=v3;return(t);}/***将矩阵值转存在一个数组里，以节省存储空间***/void storage(double C[5][501],double b,double c){ int i=0,j=0;C[i][j]=0,C[i][j+1]=0;for(j=2;j<=500;j++)C[i][j]=c;i++;j=0;C[i][j]=0;for(j=1;j<=500;j++)C[i][j]=b;i++;for(j=0;j<=500;j++)C[i][j]=(1.64-0.024*(j+1))*sin(0.2*(j+1))-0.64*exp(0.1/(j+1));i++;for(j=0;j<=499;j++)C[i][j]=b;C[i][j]=0;i++;for(j=0;j<=498;j++)C[i][j]=c;C[i][j]=0,C[i][j+1]=0;}//用于求解最大的特征值，幂法double mifa(double C[5][501]){ int m=0,i,j;double b2,b1=0,sum;double u[501],y[501];for (i=0;i<501;i++){ u[i] = 1.0;}do{ sum=0;if(m!=0)b1=b2;m++;for(i=0;i<=500;i++)sum+=u[i]*u[i];for(i=0;i<=500;i++)y[i]=u[i]/sqrt(sum);for(i=0;i<=500;i++){ u[i]=0;for(j=Max(i-r,0);j<=Min(i+s,500);j++)u[i]=u[i]+C[i-j+s][j]*y[j];}b2=0;for(i=0;i<=500;i++)b2=b2+y[i]*u[i];}while(fabs(b2-b1)/fabs(b2)>=1.0e-12);return b2;}/*****行列式LU分解*****/void LU(double C[5][501]){ double sum;int k,i,j;for(k=1;k<=501;k++){ for(j=k;j<=Min(k+s,501);j++){ sum=0;for(i=maxt(1,k-r,j-s);i<=k-1;i++)sum+=C[k-i+s][i-1]*C[i-j+s][j-1];C[k-j+s][j-1]-=sum;}for(j=k+1;j<=Min(k+r,501);j++){ sum=0;for(i=maxt(1,j-r,k-s);i<=k-1;i++)sum+=C[j-i+s][i-1]*C[i-k+s][k-1];C[j-k+s][k-1]=(C[j-k+s][k-1]-sum)/C[s][k-1];}}}/***带状DOOLITE分解，并且求解出方程组的解***/void solve(double C[5][501],double x[501],double b[501]){ int i,j,k,t;double B[5][501],c[501];for(i=0;i<=4;i++){ for(j=0;j<=500;j++)B[i][j]=C[i][j];}for(i=0;i<=500;i++)c[i]=b[i];for(k=0;k<=500;k++){ for(j=k;j<=Min(k+s,500);j++){ for(t=Max(0,Max(k-r,j-s));t<=k-1;t++)B[k-j+s][j]=B[k-j+s][j]-B[k-t+s][t]*B[t-j+s][j];}for(i=k+1;i<=Min(k+r,500);i++){ for(t=Max(0,Max(i-r,k-s));t<=k-1;t++)B[i-k+s][k]=B[i-k+s][k]-B[i-t+s][t]*B[t-k+s][k];B[i-k+s][k]=B[i-k+s][k]/B[s][k];}}for(i=1;i<=500;i++)for(t=Max(0,i-r);t<=i-1;t++)c[i]=c[i]-B[i-t+s][t]*c[t];x[500]=c[500]/B[s][500];for(i=499;i>=0;i--){ x[i]=c[i];for(t=i+1;t<=Min(i+s,500);t++)x[i]=x[i]-B[i-t+s][t]*x[t];x[i]=x[i]/B[s][i];}}//用于求解模最大的特征值，反幂法double fmifa(double C[5][501]){ int m=0,i;double b2,b1=0,sum=0,u[501],y[501];for (i=0;i<=500;i++){ [i] = 1.0;}do{ if(m!=0)b1=b2;m++;sum=0;for(i=0;i<=500;i++)sum+=u[i]*u[i];for(i=0;i<=500;i++)y[i]=u[i]/sqrt(sum);solve(C,u,y);b2=0;for(i=0;i<=500;i++)b2+=y[i]*u[i];}while(fabs(b2-b1)/fabs(b2)>=1.0e-12);return 1/b2;}/***主程序***/void main(){ double b=0.16,c=-0.064,det=1.0;int i;double C[5][501],cond;storage(C,b,c); //进行C的赋值cout.precision(12); //定义输出精度double k1=mifa(C); //利用幂法计算矩阵的最大特征值和最小特征值if(k1<0)printf("λ1=%.12e\n",k1);else if(k1>=0)printf("λ501=%.12e\n",k1);for(i=0;i<501;i++)C[2][i]=C[2][i]-k1;double k2=mifa(C)+k1;if(k2<0)printf("λ1=%.12e\n",k2);else if(k2>=0)printf("λ501=%.12e\n",k2);storage(C,b,c);double k3=fmifa(C); //利用反幂法计算矩阵A的按模最小特征值printf("λs=%.12e\n",k3);storage(C,b,c); //计算最接近特征值double u[39]={0};for(i=0;i<39;i++){ u[i]=k1+(i+1)*(k2-k1)/40;C[2][i]=C[2][i]-u[i];u[i]=fmifa(C)+u[i];printf("与数u%d 最接近的特征值λ%d: %.12e\n",i+1,i+1,u[i]);}if(k1>0) //计算矩阵A的条件数，取2范数cond=fabs(k1/k3);else if(k1<0)cond=fabs(k2/k3);storage(C,b,c);LU(C); //利用LU分解计算矩阵A的行列式for(i=0;i<501;i++)det*=C[2][i];printf("\ncond(A)=%.12e\n",cond);printf("\ndet(A)=%.12e\n",det);}三、计算结果：四、结果分析迭代初始向量的选择对果有一定的影响，选择不同的初始向量可能会得到不同阶的特征值。
数值分析作业一
习题 9：
解：求1
0(arccos )n I x dx =⎰的稳定递推公式 21/20/2/2
100n 1n 2~00001()/221101221y=x=cosy dx=-sinydy
x [0,/2]
I siny.cos .n.cosy..(/2)-n(n-1)I =(1)!E (n )(1)(n 1)!E (n n n n n n n n n y dy
y y y dy
n I E I I E n E πππππ---+∈==-+=-=-=--⎰⎰令arccosx ，则有，其中则的误差可以设为，根据误差的传递可得：
其中为偶数；同理，其中22n 2n )n 11I =I +(/2)(n 1)1
n n n π-----为奇数。
所以误差随着的变大而逐渐累加，顾不是稳定递推公式。
可求得稳定递推公式为：
习题10：
求 的稳定递推公式： 解： 1n 1011......(1)44c =-(-)44
1)c>4n 41c c 2c<4n c=n n n n n n n c I n
c E E I n
---+==+=求的递推公式为I 则有E ，根据误差的传递可得E 讨论：当时，递推公式(1)属于病态问题，即误差随增加而增加，所以递推公式要变为I ）当时，误差随增加而变小，所以递推公式（1）是稳定的
3）当4时，误差不变，递推公式（1）是稳定的。
n
n x I dx x c
1
04=+⎰
实验题
程序：
ess=input('Enter the number of ess:');
ve=zeros(1,21);
ve(2)=ess;
y=roots(poly(1:20)+ve)
plot(y)
1.当ess取值大于0.00000000001时会出现“复数”根。
表明有些解对如此扰动敏感性较大。
2.当将方程（1.2）中的扰动项改成18x 或其它形式，实验中不会出现“复数”根，各跟的抗干扰性变强。
思考题一程序如下
ess=input('Enter the number of ess:');
ve=zeros(1,21);
ve(3)=ess;
y=solve(poly2sym(poly(1:20)+ve),'x')
plot(y)
输入不同的ess值发现各根的精确度变高，干扰也变大。
思考题二程序如下：
Y=0.1;i=1;
n=input('Enter the limit value:');
while i<n
Y=Y+0.1;
i=i+1;
end
E=Y-100;
fprintf('The result of E is: %d\n',E)
思考题三程序如下：
n=1:1:10000000;
E=abs((1+n.^-1).^n-exp(1));
[Emax,nmax]=max(E)。
