(完整版)数值分析第一次作业
数值分析
数值分析第一次作业信计2 20121314044 王峥虹一、实验内容：1、已知函数在下列各点的值为：38.064.081.092.098.0|0.18.06.04.02.0|y x -------------------试用4次牛顿插值多项式)(4x P 及三次样条函数)(x S （自然边界条件）对数据进行插值，用图给出(){}10,11,1,008.02.0,=+=i i x y x i i i ，，，)(4x P 及)(x S 。
分析：先求4次插值多项式:根据差分形式的牛顿差值公式：))...(](,...,,[...))(](,,[)](,[)()(1010102100100---++--+-+=n n n x x x x x x x f x x x x x x x f x x x x f x f x Px=[0.2,0.4,0.6,0.8,1.0];y=[0.98,0.92,0.81,0.64,0.38];n=length(y);z=zeros(n,n);for i=1:nz(i,1)=y(i);endfor k=2:nfor l=k:nz(l,k)=(z(l,k-1)-z(l-1,k-1))/(x(l)-x(l-k+1));endendz结果：4次牛顿插值多项式为：)6.0)(4.0)(2.0(2083.0)4.0)(2.0(625.0)2.0(3.098.04---------=x x x x x x P )8.0)(6.0)(4.0)(2.0(5208.0-----x x x x再求三次样条插值函数：由上面及已知的：⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎣⎡---=⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎣⎡⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎣⎡075.65.475.30200005.025.00005.025.00005.025.00000243210M M M M M 程序如下：A=[2,0,0,0,0;0.5,2,0.5,0,0;0,0.5,2,0.5,0;0,0,0.5,2,0.5;0,0,0,0,2];B=[0,-3.75,-4.5,-6.75,0]';M=inv(A)*B结果：则由表达式：j j j j j j j j j j j j j j j j h x x h M y h x x h M y h x x M h x x M x S -⎪⎪⎭⎫ ⎝⎛-+-⎪⎪⎭⎫ ⎝⎛-+-+-=+++++666)(6)()(2111231311,...,1,0-=n j得，三次样条插值多项式为：⎪⎪⎩⎪⎪⎨⎧∈-+-+----∈-+-+----∈-+-+----∈-+-+---=]0.1,8.0[),8.0(9.1)0.1(3036.3)8.0(0)0.1(5893.2]8.0,6.0[),6.0(3036.3)8.0(0857.4)6.0(5893.2)8.0(8929.0]6.0,4.0[),4.0(0857.4)6.0(6536.4)4.0(8929.0)6.0(3393.1]4.0,2.0[),2.0(6536.4)4.0(9.4)2.0(3393.1)4.0(0)(3333333x x x x x x x x x x x x x x x x x x x x x S 绘制4次插值多项式及三次样条插值多项式的图像：代码：x=[0.2,0.4,0.6,0.8,1.0];y=[0.98,0.92,0.81,0.64,0.38];plot(x,y)hold onfor i=1:1:5y(i)=0.98-0.3*(x(i)-0.2)-0.625*(x(i)-0.2)*(x(i)-0.4)-0.20833*(x(i)-0.2)*(x(i)-0.4)*(x(i)-0.6)-0.5 2083*(x(i)-0.2)*(x(i)-0.4)*(x(i)-0.6)*(x(i)-0.8)endk=[0 1 10 11];x0=0.2+0.08*k;y0=zeros(4);for i=1:1:4y0(i)=0.98-0.3*(x(i)-0.2)-0.625*(x(i)-0.2)*(x(i)-0.4)-0.20833*(x(i)-0.2)*(x(i)-0.4)*(x(i)-0.6)-0. 52083*(x(i)-0.2)*(x(i)-0.4)*(x(i)-0.6)*(x(i)-0.8)endplot(x0,y0,'o',x0,y0)hold ony1=spline(x,y,x0)plot(x0,y1,'o')hold ons=csape(x,y,'variational')fnplt(x,'r')hold ongtext('原图像')gtext('三次样条自然边界')gtext('4次牛顿插值')一、实验内容：2、在区间[]11，-上分别取20,10=n 用两组等距节点对龙格函数22511)(x x f +=作多项式插值及三次样条插值，对每个n 值，分别画出插值函数及)(x f 的图形。
数值分析大作业一
数值分析大作业一一、算法设计方案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];}}}}三、程序结果。
数值分析课第一次作业答案answer1
计算机习题： 1. 作多项式 p，以 −1，0，1 为零点，首项系数为 2，并计算 p(3)。 4
答案：p = poly ([−1, 0, 1])，s = polyval(p, 3)。 2. 已知函数在下列各点的值为 xi 0.2 0.4 0.6 0.8 1.0
2
a 6 6e+154 0 1 1
b 10 10e+154 1 -1e+5 -4
c -4 -4e+154 1 1 3.999999
-1e+155 -7e+155 1e+155 答案：第二种方法更准确，因为第一种方法是一个累加的过程。 matlab 的 x = a : h : b 和 x = a + (0 : n) ∗ h 是第二种方法实现的。 代码： format long e a = 0; b = 8; n = 9; h = (b-a)/n; x(1) = a; y(1) = a; for j = 1:n, x(j+1) = x(j) + h; y(j+1) = y(1) + j*h; end [x',y',(a:h:b)',a+(0:n)’*h] 第二章 插值法 1. 当 x = 1, −1, 2 时，f (x) = 0, −3, 4，求 f (x) 的二次插值多项式。 （计算两遍，分别用拉格朗日插值和牛顿插值）
5
f (xi ) 0.98 0.92 0.81 0.64 0.38 求 4 次牛顿插值多项式 P4 (x) 并画图。 答案： 代码： x=0.2:0.2:1.0; y=[0.98,0.92,0.81,0.64,0.38]; n = length(y); if length(x)~=n, error('x and y are not compatible'); end D = zeros(n,n); D(:,1)=y(:); for j=2:n for i=j:n D(i,j) = (D(i,j-1)-D(i-1,j-1))/(x(i)-x(i-j+1)); end end p=D(1,1)*[zeros(1,n-1),1]; for k=2:n p=p+D(k,k)*[zeros(1,n-k),poly(x(1:k-1))]; end x=0.2:0.01:1.0; z=polyval(p,x); plot(x,z) 比较：p = polyf it(x, y, 4)。
JZX高等数值分析第一次实验作业
相对残差 6.1302e-16 8.7797e-09 8.0295e-09 8.5677e-09 9.1433e-09
a、 m=1 (左为相对残差，右为取对数情况)
1 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1
0 1
1.5
0
-5
-10
-15
-20
-25
-30
-35
-40
1000
4.761361
2.9675e-11
5.7069e-04
3.5336e+11
1500
15.646571
1.0778e-11
6.8236e-04
2.3672e+10
3000
132.198513
6.3164e-13
1.1700e-04
2.8110e+11
a、n=1000 时步数与相对残差关系图：（上为相对残差，下为取对数情况下结果）
（2）当 A 最大特征值远大于第二个特征值，最小特征值远小于第二个最小特征值时收敛
性情况。
思路：构造题目要求的矩阵 A。首先随机生成 n 阶矩阵 B，B 不满秩，构造对角阵 A1（最
大特征值远大于第二个最大特征值，最小特征值远小于第二个最小特征值），则由此构
造出对称正定矩阵 A： b1=B’*B; A=b1’*A1*b1。同样设定精确解 Xj 为元素全部为 1 的 n
5、 构造对称不定的矩阵，验证 Lanczos 方法的近似中断，观察收敛曲线中的峰点个数和特
征值的分布关系；观察当出现峰点时，MINRES 方法的收敛性态怎样。
解：思路：类似前两题，首先构造出一个 n 阶对角阵 D，其对角线上有 m 个负值，再对随
数值分析第一次作业
《数值分析》计算作业院系：航空科学与工程学院学号： SY1005512姓名：王天龙日期： 2010年10月31日计算实习说明书目的：训练运用计算机进行科学与工程计算的能力。
要求：1．独立进行算法设计、程序设计和上机运算，并得出正确的结果。
2．编制程序时全部采用双精度，要求按题目的要求设计输出，并执行打印。
3．只能根据题目给出的信息并且只允许一次计算得出全部结果。
题目：第一题 设有501×501的矩阵123499500501a b c b a b cc b a b c A c b a b c c b a b c ba ⎡⎤⎢⎥⎢⎥⎢⎥⎢⎥=⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎣⎦其中0.1(1.640.024)s i n (0.2)0.64 (125i i a i i e i =--= ，，，；0.16b =；0.064c =-。
矩阵A 的特征值12501()ii λ= ，，，满足 125011501||min ||S i i λλλλλ≤≤<<<= ，试求：1.1λ，501λ和S λ的值。
2.A 的与数5011140k kλλμλ-=+最接近的特征值(1239)ik k λ= ，，，。
3.A 的（谱范数）条件数2()cond A 和行列式det A 。
说明：1．在所有的算法中，凡是要给出精度水平ε的，都取1210ε-=。
2．选择算法时，应使A 的所有零元素都不存储。
3．打印以下内容： （1）算法的设计方案。
（2）全部源程序（要求注明主程序和每个子程序的功能）。
（3）特征值1λ，501λ，S λ和(1239)ik k λ= ，，，以及2()cond A ，det A 的值。
4.采用e 型输出所有计算结果，并至少显示12位有效数字。
一、程序算法的设计算法设计方案如下：二、全部源程序编程软件：Fortran：三、计算结果1.特征值1λ，501λ，S λ1-.107001135582E+02λ=，501 .972463398616E+01λ=，-.555823879237E-02S λ=2. (1239)ik k λ= ，，，如下表所示（ZK 代表ik λ）3.A 的条件数2()cond A 和行列式det A 的值2() .192509100000E+04cond A =，det .277059968428+119A =五、讨论这里选取的初始向量为X(i)=1,X={x1,x2,x3,,,,,,,x501},当初始向量与特征向量较近时，收敛较快，若初始向量与特征向量正交，则求解可能失真。
清华大学高等数值分析 第一次实验作业
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
数值分析第一次作业解答
数值分析第一次作业解答1：(a) —个问题的病态性如何，与求解它的算法有关系。
x ;(b) 无论问题是否病态，好的算法都会得到它好的近似解。
x ；(C)计算中使用更高的精度，可以改善问题的病态性。
X ；(d) 用一个稳定的算法计算一个良态问题，一定会得到他好的近似解。
V;(e) 浮点数在整个数轴上是均匀分布。
x ;(f) 浮点数的加法满足结合律。
x(g) 浮点数的加法满足交换律。
X ;(h) 浮点数构成有效集合。
V;(i) 用一个收敛的算法计算一个良态问题，一定得到它好的近似解2: 程序t=0.1;n=1:10;e=n/10-n*te = 1.0e-015 *[ 0 0 -0.0555 0 0-0.1110 -0.1110 0 0 0] 由舍人误差造成n=3，6，7 时的结果不为零。
4:两种等价的一元二次方程求解公式-b - Pb2 - 4acx =2a2cx 二-b b2 - 4ac对a=1, b=-100000000, c=1，应采用哪种算法?A二[1,-100000000,1];roots(A);可得：X1 = 100000000;x2=0a=1;b=-100000000;c=1;x1仁(-b-sqrt(b*b-4*a*c))/(2*a)x12=(-b+sqrt(b*b-4*a*c))/(2*a)x2仁2*c/(-b-sqrt(b*b-4*a*c))x22=2*c/(-b+sqrt(b*b-4*a*c))由第一种算法：X1 = 100000000;x2=7.45058 X10由第二种算法：X1 = 13417728;x2=-1.0 X108原因：太小的数作分母。
5:程序:fun cti on y=tt(x)s=0;t=x;n=1;while s+t~=s;s=s+t;t=-x A2/(( n+1)*( n+2))*tn=n+2;endntt(2n 1)eps)(a)t小于计算机的计算精度。
数值分析作业题(1)
第一章 误差与算法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：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. 数值计算方法设计的基本手段是（）.（A ）近似 （B ） 插值 （C ） 拟合 （D ） 迭代2. 为了在有限时间内得到结果，用有限过程取代无限过程所产生的近似解与精 确解之间的误差称为（）.（A ）舍入误差 （B ） 截断误差 （C ） 测量误差 （D ） 绝对误差3. 由于计算机的字长有限，原始数据在机器内的表示以及进行算术运算所产生 的误差统称为（）.（A ）舍入误差 （B ） 截断误差 （C ）相对误差 （D ） 绝对误差 4. 数值计算方法研究的核心问题可以概括为（）对计算结果的影响. （A ）算法的稳定性（B ） 算法的收敛性 （C ）算法的复杂性 （D ） 近似 N dx5. 当N 充分大时,利用下列各式计算I 二 半,等式（）得到的结果最好. •N 1 +x（A ） I =arcta n （N 1）-arcta n （N ） （B ）I 二 arcta n （N 2 N 1） 6.计算（、、2-1）6,取；2 1.4,利用下列哪个公式得到的结果最好 ？为什么？7. 计算球体的体积，已知半径的相对误差限不超过 3 10”，则计算所得体积的相 对误差限如何估计？8. 设x 0,近似值x 的相对误差限为:,试估计In x 的误差限.9. 计算圆柱体的体积，已知底面半径r 及圆柱高h 的相对误差限均不超过、：，则 计算所得体积的相对误差限如何估计?.10. 用秦九韶算法求f （x ） = 4x 3「3x 2 • x 「1在x = 2处的值.111. 已知近似值 X =1.0000 的误差限；（x ）=1 10，, f （x ）二丄 X 2,求；（f （x ））,并16说明X”及f （X”）的各有几位有效数字.(C)I = arcta n( 2 N 2 N 1) (D)(A) 1(；2 1)6 (B) (3-2运3 (C) 1(3 - 2 ・(D) 99 - 70、212. 设a为非零常数，已知y0的近似值y0,由递推式y n =ay n斗计算序列｛y n｝的近似值,分析该算法的稳定性.。
数值分析第一次大作业
《数值分析》计算实习报告第一题院系：机械工程及自动化学院_学号： _____姓名： _ ______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=时，计算的结果如下：此结果即为上文中的正确计算结果。
数值分析第一次作业答案
作业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、最大特征值λ501，最小特征值λ1，按模最小特征值λs的计算方法首先使用一次幂法运算可以得到矩阵的按模最大的特征值λ，λ必为矩阵A的最大或最小特征值，先不做判断。
对原矩阵A进行一次移项，即（A-λI），在进行一次幂法运算，可以得到另一个按模最大特征值λ0。
比较λ和λ的大小，较大的即为λ501，较小的即为λ1。
对矩阵A进行一次反幂法运算，即可得到按模最小特征值λs。
2、A与μk 值最接近的特征值λik的计算方法首先计算出k所对应的μk 值，对原矩阵A进行一次移项，即（A-μkI），得到一个新的矩阵，对新矩阵进行一次反幂法运算，即可得到一个按模最小特征值λi 。
则原矩阵A与μk值最接近的特征值λik=λi+μk。
3、A的（谱范数）条件数cond（A）2的计算方法其中错误!未找到引用源。
矩阵A的按模最大和按模最小特征值。
（二）程序编写思路。
由于算法要求A的零元素不存储，矩阵A本身为带状矩阵，所以本题的赋值，LU分解，反幂法运算过程中，均应采用Doolittle分解法求解带状线性方程组的算法思路。
幂法、反幂法和LU分解均是多次使用，应编写子程序进行反复调用。
二、源程序：#include<stdio.h>#include<iostream>#include<stdlib.h>#include<math.h>#include<float.h>#include<iomanip> /*头文件*//*定义全局变量*/#define N 502 /*取N为502，可实现从1到501的存储，省去角标变换的麻烦*/ #define epsilon 1.0e-12 /*定义精度*/#define r 2 /*r，s为带状矩阵的半带宽，本题所给矩阵二者都是2*/ #define s 2double c[6][N]; /*定义矩阵存储压缩后的带状矩阵*/double fuzhi(); /*赋值函数*/void LUfenjie(); /*LU分解程序*/int max(int a,int b); /*求两个数字中较大值*/int min(int a,int b); /*求两个数字中较小值*/double mifa(); /*幂法计算程序*/double fanmifa(); /*反幂法计算程序*/double fuzhi() /*赋值程序，按行赋值，行从1到5，列从1到501，存储空间的第一行第一列不使用，角标可以与矩阵一一对应，方便书写程序*/{int i,j;i=1;for(j=3;j<N;j++){c[i][j]=-0.064;}i=2;for(j=2;j<N;j++){c[i][j]=0.16;}i=3;for(j=1;j<N;j++){c[i][j]=(1.64-0.024*j)*sin(0.2*j)-0.64*exp(0.1/j);}i=4;for(j=1;j<N-1;j++){c[i][j]=0.16;}i=5;for(j=1;j<N-2;j++){c[i][j]=-0.064;}return(c[i][j]);}int max(int a,int b){ return((a>b)?a:b);}int min(int a,int b){ return((a<b)?a:b);}void LUfenjie() /*LU分解程序，采用的是带状矩阵压缩存储后的LU分解法*/{double temp;int i,j,k,t;for(k=1;k<N;k++){ for(j=k;j<=min(k+s,N-1);j++){temp=0;for(t=max(1,max(k-r,j-s));t<=(k-1);t++){temp=temp+c[k-t+s+1][t]*c[t-j+s+1][j];}c[k-j+s+1][j]=c[k-j+s+1][j]-temp;}for(i=k+1;i<=min(k+r,N-1);i++){temp=0;for(t=max(1,max(i-r,k-s));t<=(k-1);t++){temp=temp+c[i-t+s+1][t]*c[t-k+s+1][k];}c[i-k+s+1][k]=(c[i-k+s+1][k]-temp)/c[s+1][k];}}}double mifa() /*幂法计算程序*/ {double u0[N],u1[N];double temp,Lu,beta=0,beta0;int i,j;for(i=1;i<N;i++) /*选取迭代初始向量*/{u0[i]=1;}do{beta0=beta;temp=0;for(i=1;i<N;i++){temp=temp+u0[i]*u0[i]; }Lu=sqrt(temp);for(i=1;i<N;i++){u1[i]=u0[i]/Lu;}for(i=1;i<N;i++){temp=0;for(j=max(i-1,1);j<=min(i+2,N-1);j++){temp=temp+c[i-j+s+1][j]*u1[j]; }u0[i]=temp;} //新的u0temp=0;for(i=1;i<N;i++){temp=temp+u1[i]*u0[i]; }beta=temp;}while(fabs(beta-beta0)/fabs(beta)>=epsilon); /*迭代运行条件判断*/return(beta);}double fanmifa() /*反幂法计算程序*/{double u0[N],u1[N],u2[N];double temp,Lu,beta=0,beta0;int i,t;for(i=1;i<N;i++) /*选取迭代初始向量*/{u0[i]=1;}do{beta0=beta;temp=0;for(i=1;i<N;i++){temp=temp+u0[i]*u0[i]; }Lu=sqrt(temp);for(i=1;i<N;i++){u1[i]=u0[i]/Lu;u2[i]=u1[i];}fuzhi();LUfenjie();/*带状矩阵压缩存储并进行LU分解后，求解线性方程组得到迭代向量u k，即程序中的u0*/for(i=2;i<N;i++){ temp=0;for(t=max(1,i-r);t<=(i-1);t++){temp=temp+c[i-t+s+1][t]*u2[t];}u2[i]=u2[i]-temp;}u0[N-1]=u2[N-1]/c[s+1][N-1];for(i=N-2;i>=1;i--){ temp=0;for(t=i+1;t<=min(i+s,N-1);t++){temp=temp+c[i-t+s+1][t]*u0[t];}u0[i]=(u2[i]-temp)/c[s+1][i];}temp=0;for(i=1;i<N;i++){temp=temp+u1[i]*u0[i]; }beta=temp;beta=1/beta; /*beta即为所求特征值，可直接返回*/}while(fabs(beta-beta0)/fabs(beta)>=epsilon); /*迭代运行条件判断*/return(beta);}void main(){double u[40]; /*定义数组，存放k值运算得到的μk值*/double lambda1,lambda501,lambdak,a,b,d,cond,det;int i,j,k;fuzhi();a=mifa(); /*幂法计算按模最大值*/fuzhi();d=fanmifa(); /*反幂法计算按模最小值*/fuzhi();for(j=1;j<N;j++){c[3][j]=c[3][j]-a;}b=mifa()+a; /*移项后幂法计算按模最大值*/if(a>b) /*比较两个按模最大值大小，并相应输出最大特征值λ501和最小特征值λ1*/ {lambda1=b;lambda501=a;printf("矩阵A最小的特征值lambda1=%13.11e\n",lambda1);printf("矩阵A最大的特征值lambda501=%13.11e\n",lambda501);}else{lambda1=a;lambda501=b;printf("矩阵A最小的特征值lambda1=%13.11e\n",lambda1);printf("矩阵A最大的特征值lambda501=%13.11e\n",lambda501);}printf("矩阵A按模最小特征值lambdas=%13.11e\n",d); /*输出按模最小特征值λs*/for(k=1;k<40;k++) /*对每一个进行移项反幂法运算，求出最接近μk的特征值并输出*/ {u[k]=(lambda501-lambda1)*k/40+lambda1;fuzhi();for(j=1;j<N;j++){c[3][j]=c[3][j]-u[k];}lambdak=fanmifa()+u[k];i=k;printf("矩阵A最接近uk特征值lambdak%d=%13.11e\n",i,lambdak);}cond=fabs(a/d);printf("A的条件数=%13.11e\n",cond); /*计算A条件数并输出*/fuzhi(); /*计算A的行列式值并输出*/LUfenjie();det=1;for(i=1;i<N;i++){det=det*c[3][i];}printf("行列式的值detA=%13.11e\n",det);}三、程序的运行结果：四、初始向量的选取对计算结果的影响：（一）选取形式不变，数值变换1、取u0为[0.5,0.5………..0.5]，运行结果如下：2、取u0为[50,50………..50]，运行结果如下：从运行结果来看，此类初始向量的选取对结果不会产生影响，即使选成0，结果也不变化。
数值分析第一次作业
数值分析第一次作业班级学号姓名习题24、用Newton法求方程f(x)=x^3-2*x^2-4*x-7=0在[3，4]中的根。
代码：function[x_star,k]=Newton1[fname,dfname,x0,ep,Nmax]if nargin<5 Nmax=500; endif nargin<4 ep=1e-5;endx=x0;x0=x+2*ep;k=0;while abs(x0-x)>ep&k<Nmax k=k+1x0=x;x=x0-feval(fname,x0)/feval(dfname,x0);endx_star=x;if k==Nmax warning(‘已迭代上限次数’);endfname=inline('x^3-2*x^2-4*x-7');dfname=inline('3*x^2-4*x-4');[x_star,k]=Newton1(fname,dfname,3.5)x_star =3.6320k =4方法二：2－4用割线法求方程的根function [x_star,k]=Gline(fun,x0,x1,ep,Nmax)if nargin<5 Nmax=500;endif nargin<4 ep=1e-5;endk=0;while abs(x1-x0)>ep&k<Nmaxk=k+1;x2=x1-feval(fun,x1)*(x1-x0)/(feval(fun,x1)-feval(fun,x0))x0=x1;x1=x2;endx_star=x1;if k==Nmax warning('已迭代上限次数');endfun=inline('x^3-2*x^2-4*x-7');[x_star,k]=Gline(fun,3,4)x2 =3.5263x2 =3.6168x2 =3.6327x2 =3.6320x2 =3.6320x_star =3.6320k =5习题33、用列主元消去法解方程组[-1 2 -2; 3 -1 4; 2 -3 -2][x1 x2 x3]=[-1 7 0] 代码：function x=Gauss_x1(A,b)A=[A’;b]’,n=length(b);for k=1:n-1s=A(k,k);p=k;for i=l+1:nif abs(s)<abs(A(i,k))s=A(I,k);p=I;endendAfor i=k+1:nm=A(i,k)/A(k,k);fprintf(‘m%d%d=%f\n’,i,k,m);for j=k:n+1A(i,j)=A(i,j)-m*A(k,j);endendfprintf(‘A%d=\n’,k+1);AendA(n,n+1)=A(n,n+1)/A(n.n);for i=n-1:-1:1s=0for j=i+1:ns=s+A(i,j)*A(j,n+1);endA(i,n+1)=(A(i,n+1)-s)/A(i,i);endA(:,n+1)A=[-1,2,-2;3,-1,4;2,-3,-2];b=[-1;7;0];x=Gauss_x1(A,b)A =3.0000 -1.00004.0000 7.00000 1.6667 -0.6667 1.33330 -2.3333 -4.6667 -4.6667A=3.0000 -1.00004.0000 7.00000 -2.3333 -4.6667 -4.66670 0 -4.0000 -2.0000x =2.00001.00000.50004、用追赶法解三对角方程[2 -1 0 0 0;-1 2 -1 0 0;0 -1 2 -1 0;0 0 -1 2 -1;0 0 0 -1 2][x1 x2 x3 x4 x5=[1 0 0 0 0 ] 代码：function x=zhuigan(A,B,C,D)n=length(B);Xzeros(1,n);U=zeros(1,n);Q=zeros(1,n);U(1)=C(1)/B(1);Q(1)=D(1)/B(1);for i=2:n-1U(i)=C(i)/(B(i)-U(i-1)*A(i-1));endfor i=2:nQ(i)=(D(i)-Q(i-1)*A(i-1))/(B(i)-U(i-1)*A(i-1));endX(n)=Q(n);for i=n-1:-1:1X(i)=Q(i)-U(i)*X(i+1);endXA=[-1，-1，-1，-1;];B=[2,2,2,2,2];C=[-1,-1,-1,-1];D=[1;0;0;0;0];X=zhuigan(A,B,C,D)X= 0.8333 0.6667 0.5000 0.3333 0.16676、用三角分解法解方程组[-2 4 8;-4 18 -16;-6 2 -20][x1 x2 x3]=[5 8 7]代码function[y,x]=LU_s(A,b)b=b';A=[A';b]',n=length(b');x=zeros(n,1);y=zeros(n,1);U=zeros(n);L=eye(n);for k=1:nU(1,k)=A(1,k);L(k,1)=A(k,1)/U(1,1);endfor i=2:nfor k=i:nlu=0;lu1=0;for j=1:i-1lu=lu+L(i,j)*U(j,k);lu1=lu1+L(k,j)*U(j,i);endU(i,k)=A(i,k)-lu;L(k,i)=(A(k,i)-lu1)/U(i,i);endendLUfor i=1:nly=0;for j=1:ily=ly+L(i,j)*y(j);endy(i)=b(i)-ly;endfor i=n:-1:1ly1=0;for j=i+1:nly1=ly1+U(i,j)*x(j);endx(i)=(y(i)-ly1)/U(i,i);endA=[-2,4,8;-4,18,-16;-6,2,-20];b=[5;8;7];[y,x]=LU_s(A,b)A =-2 4 8 5-4 18 -16 8-6 2 -20 7L =1 0 02 1 03 -1 1U =-2 4 80 10 -320 0 -76y =5-2-10x =-1.53160.22110.13163-8用LU分解法解线性方程组[5,7,9,10;6,8,10,9;7,10,8,7;5,7,6,5][x1 x2 x3 x4]=[1 1 1 1] 代码function[y,x]=LU_s(A,b)b=b';A=[A';b]',n=length(b');x=zeros(n,1);y=zeros(n,1);U=zeros(n);L=eye(n);for k=1:nU(1,k)=A(1,k);L(k,1)=A(k,1)/U(1,1);endfor i=2:nfor k=i:nlu=0;lu1=0;for j=1:i-1lu=lu+L(i,j)*U(j,k);lu1=lu1+L(k,j)*U(j,i);endU(i,k)=A(i,k)-lu;L(k,i)=(A(k,i)-lu1)/U(i,i);EndendLUfor i=1:nly=0;for j=1:ily=ly+L(i,j)*y(j);endy(i)=b(i)-ly;endfor i=n:-1:1ly1=0;for j=i+1:nly1=ly1+U(i,j)*x(j);endx(i)=(y(i)-ly1)/U(i,i);endA=[5,7,9,10;6,8,10,9;7,10,8,7;5,7,6,5];b=[1;1;1;1];[y,x]=LU_s(A,b)A =5 7 9 10 16 8 10 9 17 10 8 7 15 76 5 1L =1.0000 0 0 01.2000 1.0000 0 01.4000 -0.5000 1.0000 01.0000 0 0.6000 1.0000U =5.0000 7.0000 9.0000 10.0000 0 -0.4000 -0.8000 -3.0000 0 0 -5.0000 -8.5000 0 0 0 0.1000y =1.0000 -0.2000 -0.5000 0.3000x =20.0000 -12.0000 -5.0000 3.0000。
北航数值分析第一次大作业(幂法反幂法)
一、问题分析及算法描述1. 问题的提出：（1）用幂法、反幂法求矩阵A =[a ij ]20×20的按摸最大和最小特征值，并求出相应的特征向量。
其中 a ij ={sin (0.5i +0.2j ) i ≠j 1.5cos (i +1.2j ) i =j要求：迭代精度达到10−12。
（2）用带双步位移的QR 法求上述的全部特征值，并求出每一个实特征值相应的特征向量。
2. 算法的描述：(1) 幂法幂法主要用于计算矩阵的按摸为最大的特征值和相应的特征向量。
其迭代格式为：{ 任取非零向量u 0=(h 1(0)，⋯，h n (0))T|h r (k−1)|=max 1≤j≤n |h r (k−1)| y ⃑ k−1=u ⃑ k−1|h r (k−1)| u ⃑ k =Ay ⃑ k−1=(h 1(k )，⋯，h n (k ))T βk =sgn (h r (k−1))h r (k ) (k =1，2，⋯) 终止迭代的控制选用≤ε。
幂法的使用条件为n ×n 实矩阵A 具有n 个线性无关的特征向量x 1，x 2，⋯，x n ，其相应的特征值λ1，λ2，⋯，λn 满足不等式|λ1|>|λ2|≥|λ3|≥⋯≥|λn |或λ1=λ2=⋯=λm|λ1|>|λm+1|≥|λm+2|≥⋯≥|λn |幂法收敛速度与比值|λ2λ1|或|λm+1λ1|有关，比值越小，收敛速度越快。
(2) 反幂法反幂法用于计算n ×n 实矩阵A 按摸最小的特征值，其迭代格式为：{任取非零向量u 0∈R nηk−1=√u ⃑ k−1T u ⃑ k−1 y ⃑ k−1=u ⃑ k−1ηk−1⁄ Au ⃑ k =y ⃑ k−1 βk =y ⃑ k−1u ⃑ k (k =1，2，⋯) 每迭代一次都要求解一次线性方程组Au ⃑ k =y ⃑ k−1。
当k 足够大时，λn ≈1βk ，y ⃑ k−1可近似的作为矩阵A 的属于λn 的特征向量。
最新(完美版)数值分析上机作业
一. 上机作业任务一： 用五点差分格式求解Poisson 方程边值问题，P124-------3、4（任选一题）。
二. 上机作业任务二：用Simpson 求积法计算定积分()baf x dx ⎰。
下面两种方法任选一。
（1）变步长复化Simpson 求积法。
（2）自适应Simpson 求积法三. 上机作业任务三：用MATLAB 语言编写连续函数最佳平方逼近的算法程序（函数式M 文件）。
要求算法程序可以适应不同的具体函数，具有一定的通用性。
并用此程序进行数值试验，写出计算实习报告。
所编程序具有以下功能：1. 用Lengendre 多项式做基，并适合于构造任意次数的最佳平方逼近多项式。
可利用递推关系 0112()1,()()(21)()(1)()/2,3,.....n n n P x P x xP x n xP x n P x n n --===---⎡⎤⎣⎦=2. 被逼近函数f(x)用M 文件建立数学函数。
这样，此程序可通过修改建立数学函数的M 文件以适用不同的被逼近函数（要学会用函数句柄）。
3. 要考虑一般的情况]1,1[],[)(+-≠∈b a x f 。
因此，程序中要有变量代换的功能。
4. 计算组合系数时，计算函数的积分采用数值积分的方法。
5. 程序中应包括帮助文本和必要的注释语句。
另外，程序中也要有必要的反馈信息。
6. 程序输入：（1）待求的被逼近函数值的数据点0x (可以是一个数值或向量)（2）区间端点：a,b 。
7. 程序输出：（1）拟合系数：012,,,...,n c c c c（2）待求的被逼近函数值00001102200()()()()()n n s x c P x c P x c P x c P x =++++ 8． 试验函数：()cos ,[0,4]f x x x x =∈+；也可自选其它的试验函数。
9． 用所编程序直接进行计算，检测程序的正确性，并理解算法。
10． 分别求二次、三次、。
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
问题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代表端点的一阶导。
x=[0.25 0.30 0.39 0.45 0.53];y=[0.5 0.5477 0.6245 0.6708 0.7280];n=5,s=1.0,t=0.6868for j=1:1:n-1h(j)=x(j+1)-x(j);endfor j=2:1:n-1r(j)=h(j)/(h(j)+h(j-1));endfor j=1:1:n-1u(j)=1-r(j);endfor j=1:1:n-1f(j)=(y(j+1)-y(j))/h(j);endfor j=2:1:n-1d(j)=6*(f(j)-f(j-1))/(h(j-1)+h(j));endd(1)=6*(f(1)-s)/h(1)d(n)=6*(t-f(n-1))/h(n-1)a=zeros(n,n);for j=1:1:na(j,j)=2;endr(1)=1;u(n)=1;for j=1:1:n-1a(j+1,j)=u(j+1);a(j,j+1)=r(j);endb=inv(a)m=b*d'p=zeros(n-1,4); %p矩阵为S(x)函数的系数矩阵for j=1:1:n-1p(j,1)=m(j)/(6*h(j));p(j,2)=m(j+1)/(6*h(j));p(j,3)=(y(j)-m(j)*(h(j)^2/6))/h(j);p(j,4)=(y(j+1)-m(j+1)*(h(j)^2/6))/h(j);end对于第二中边界，已知边界二阶倒数，可写出下面的矩阵：⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎣⎡=⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎣⎡⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎣⎡432104321034322110d M M M M M 2000200002002002d 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 =λ0=0 d 0=d n =0解：由matlab 计算得：由此得矩阵形式的线性方程组为：⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎣⎡=⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎣⎡⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎣⎡0 2.4286-3.2667-4.3143-0M M M M M 25714.00000204286.000004000.026000.0006429.023571.000243210解得M 0=0 ;M 1= -1.8795;M 2= -0.8636; M 3= -1.0292; M 4=0S(x)=⎪⎪⎩⎪⎪⎨⎧∈-+-+--∈-+-+---∈-+-+---∈-+-+--]53.0,45.0[x )45.0(1000.953.03987.845.00-x .5302.1442-]45.0,39.0[x 9.30x 1903.11x 54.0417.1093.0-x 8590.2x .4503990.2]39.0,30.0[x 03.0x 9518.6x 9.301137.603.0-x 5993.1x .3904806.3]30.0,25.0[x )25.0(9697.1030.0 010.25-x 2652.6x 0.30 033333333，）（）（）（），（）（）（）（），（）（）（）（，）（）（）（x x x x x matlab 程序代码如下：function tgsanci(n,s,t) %n 代表元素数， x=[0.25 0.30 0.39 0.45 0.53];y=[0.5 0.5477 0.6245 0.6708 0.7280];n=5for j=1:1:n-1h(j)=x(j+1)-x(j);endfor j=2:1:n-1r(j)=h(j)/(h(j)+h(j-1));endfor j=1:1:n-1u(j)=1-r(j);endfor j=1:1:n-1f(j)=(y(j+1)-y(j))/h(j);endfor j=2:1:n-1d(j)=6*(f(j)-f(j-1))/(h(j-1)+h(j));endd(1)=0d(n)=0a=zeros(n,n);for j=1:1:na(j,j)=2;endr(1)=0;u(n)=0;for j=1:1:n-1a(j+1,j)=u(j+1);a(j,j+1)=r(j);endb=inv(a)k=am=b*d'p=zeros(n-1,4); %p矩阵为S(x)函数的系数矩阵for j=1:1:n-1p(j,1)=m(j)/(6*h(j));p(j,2)=m(j+1)/(6*h(j));p(j,3)=(y(j)-m(j)*(h(j)^2/6))/h(j);p(j,4)=(y(j+1)-m(j+1)*(h(j)^2/6))/h(j);endp由下面的一段程序，画出自然边界与第一边界的图形：程序代码如下：x=[0.25 0.30 0.39 0.45 0.53];y=[0.5 0.5477 0.6245 0.6708 0.7280];s=csape(x,y,'variational')fnplt(s,'r')hold onxlabel('x')ylabel('y')gtext('三次样条自然边界')plot(x,y,'o',x,y,':g')hold ons.coefsr=csape(x,y,'complete',[1.0 0.6868])fnplt(r,'b')gtext('三次样条第一边界')hold on r.coefs图中的三条线几乎重合。
xy图20-1 在一小段区间内的图形xy图20-2 在整个给出的区间的图形问题2：1已知函数在下列各点的值为试用4次牛顿插值多项式P 4（x ）及三次样条函数S （x ）（自然边界条件）对数据进行插值。
用图给出{（x i ，y i ），x i =0.2+0.08i ，i=0，1, 11, 10}，P 4（x ）及S （x ）。
分析：（1）要用4次牛顿插值多项式处理数据， P n =f(x 0)+f[x 0,x 1](x-x 0)+ f[x 0,x 1,x 2](x-x 0) (x-x 1)+···+ f[x 0,x 1，···x n ](x-x 0) ···(x-x n-1) 首先我们要知道牛顿插值多项式的系数，即均差表中得部分均差。
解：由所以有四次插值牛顿多项式为：P 4（x ）=0.98-0.3(x-0.2)-0.62500 (x-0.2)(x-0.4) -0.20833 (x-0.2)(x-0.4)(x-0.6)-0.52083(x-0.2)(x-0.4)(x-0.6)（x-0.8）计算牛顿插值中多项式系数的程序如下： function varargout=newtonliu(varargin) clear,clcx=[0.2 0.4 0.6 0.8 1.0];fx=[0.98 0.92 0.81 0.64 0.38]; newtonchzh(x,fx);function newtonchzh(x,fx) %由此函数可得差分表 n=length(x);fprintf('*****************差分表*****************************\n'); FF=ones(n,n); FF(:,1)=fx'; for i=2:nfor j=i:nFF(j,i)=(FF(j,i-1)-FF(j-1,i-1))/(x(j)-x(j-i+1)); end endfor i=1:nfprintf('%4.2f',x(i)); for j=1:ifprintf('%10.5f',FF(i,j)); endfprintf('\n'); end（2）用三次样条插值函数由上题分析知,要求各点的M 值： 有matlab 编程计算求出个三次样条函数：⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎣⎡=⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎣⎡⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎣⎡06.7500-4.5000-3.7500-0M M M M M 2.5000000020.50000000.50002.500000.50002.5000000 243210 解得：M 0=0 ;M 1= -1.6071;M 2= -1.0714; M 3= -3.1071; M 4=0⎪⎪⎩⎪⎪⎨⎧∈-+-+--∈-+-+-∈-+-+--∈-+-+---]0.1,8.0[x )8.0(9.10.1 3036.38.00-x 0.12.5893-]8.0,6.0[x 6.0x 3036.3x 8.00857.4.60-x 5893.2-x .800.8929-]6.0,4.0[x 4.0x 7508.4x 6.04.6536 .40-x 8929.0x .601.3393- ]4.0,2.0[x )2.0(6536.44.0900.42.0 1.33934.00 3333333，）（）（）（），（）（）（）（），（）（）（）（，）（）（）（x x x x x x x 三次样条插值函数计算的程序如下：function tgsanci(n,s,t) %n 代表元素数，s,t 代表端点的一阶导。